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Chapter  1 

Introduction  and  Summary 


This  part  of  the  ONR-funded  project  was  principally  devoted  to  the  development  of  a  simplified  time-domain 
nonlinear  model  of  thermoacoustic  devices. 

The  starting  point  for  the  development  of  this  model  is  the  observation  that,  in  most  cases,  the  size  of 
these  devices  in  the  direction  normal  to  the  wave  fronts  is  usually  greater  than  the  length  scales  parallel 
to  the  wave  fronts.  This  circumstance  suggests  the  use  of  a  well-known  approximation  in  Fluid  Mechanics, 
consisting  in  the  averaging  of  the  conservation  equations  over  the  cross  section. 

A  difficulty  inherent  in  this  procedure  is  that  averaging  removes  information  of  the  exchange  processes  of 
momentum  and  energy  between  the  working  fluid  and  the  solid  structure  in  contact  with  it.  These  processes 
are  of  crucial  importance  in  thermoacoustic  phenomena  and,  therefore,  it  is  essential  to  reintroduce  the  lost 
information  into  the  mathematical  model,  albeit  approximately. 

This  has  been  the  greatest  difficulty  encountered  in  this  work,  but  we  believe  to  have  found  a  satisfactory 
solution  as  described  in  Chapter  2.  This  chapter  contains  a  detailed  exposition  of  the  model  and  its  rationale 
and  presents  several  numerical  examples  which  demonstrate  its  versatility:  a  prime  mover,  an  externally 
driven  thermoacoustic  refrigerator,  and  a  combined  prime  mover/refrigerator  system.  In  addition  to  the 
analytical  formulation  of  the  model,  we  have  also  developed  a  robust  numerical  method  for  the  soluton  of 
the  equations. 

The  outcome  of  this  work  is  therefore  a  complete  (if  approximate)  fully  nonlinear  time-domain  model  of 
thermoacoustic  devices  -  the  first  of  its  kind.  While  certainly  incomplete,  the  model  agrees  reasonably  well 
with  data,  is  numerically  efficient  and,  at  the  very  least,  it  is  useful  to  compare  design  options. 

In  parallel  with  this  work,  we  have  also  developed  a  time-domain  formulation  of  the  linear  theory  of 
thermoacoustic  devices,  as  well  as  a  weakly  non-linear  theory;  these  results  are  described  in  chapter  3. 

The  last  part  of  the  study  concerns  a  different  topic  -  the  propagation  of  pressure  waves  through  layers  of 
bubbly  liquids.  The  motivcation  of  the  work  is  the  possible  use  of  such  systems  as  parametric  arrays  for  the 
generation  of  low-frequency  sound,  and  also  for  the  attenuation  of  strong  waves;  these  results  are  described 
in  chapter  4. 

1.1  Publications 

The  work  has  given  rise  to  the  following  publications,  which  are  here  listed  with  the  respective  abstracts: 

Watanabe,  M.,  Prosperetti,  A.  &  Yuan,  H.,  A  Simplified  Model  for  Linear  and  Nonlinear 
Processes  in  Thermoacoustic  Prime  Movers.  Part  I.  Model  and  linear  theory.  J.  Acoust.  Soc . 
Am.  102,  3484-3496,  1997 

ABSTRACT:  A  simplified  quasi-one-dimensional  model  of  thermoacoustic  devices  is  formulated  by  averaging 
the  conservation  equations  over  the  cross  section.  Heat  transfer  and  drag  effects  are  introduced  by  means  of 
suitable  coefficients.  While  the  primary  motivation  for  this  work  is  the  development  of  a  model  approximately 


valid  in  the  nonlinear  regime,  the  focus  of  this  paper  is  the  proper  formulation  of  the  transfer  coefficients  and 
the  analysis  of  the  linear  problem.  The  accuracy  of  the  model  is  demonstrated  by  comparison  with  existing 
more  precise  theories  and  data.  Examples  of  devices  with  variable  cross  section  demonstrate  the  flexibility 
of  the  approach. 

Yuan,  H.,  Karpov,  S.  &  Prosperetti,  A.,  A  Simplified  Model  for  Linear  and  Nonlinear  Processes 
in  Thermoacoustic  Prime  Movers.  Part  II.  Nonlinear  Oscillations.  J.  Acoust.  Soc.  Am,  102, 
3497-3506,  1997 

ABSTRACT:  The  simplified  quasi-one-dimensional  model  of  thermoacoustic  devices  formulated  in  Part  I 
(Watanabe,  Prosperetti,  &  Yuan,  J.  Acoust.  Soc.  Am.  102,  3484-3496,  1997)  is  studied  in  the  nonlinear 
regime.  A  suitable  numerical  method  is  described  which  is  able  to  deal  with  the  steep  waveforms  that  develop 
in  the  system  without  inducing  spurious  oscillations,  appreciable  numerical  damping,  or  numerical  diffusion. 
The  results  are  compared  with  some  experimental  ones  available  in  the  literature.  Several  of  the  observed 
phenomena  are  reproduced  by  the  model.  Quantitative  agreement  is  also  reasonable  when  allowance  is  made 
for  likely  temperature  non-uniformities  across  the  heat  exchangers. 

Karpov,  S.  &:  Prosperetti,  A.,  Linear  Thermoacoustic  Instability  in  the  Time  Domain.  J. 
Acoust  Soc.  Am.  103,  3309-3317,  1998 

ABSTRACT:  An  approximate  time-domain  analysis  of  the  development  of  the  thermoacoustic  instability  in 
gas-filled  tubes  is  developed  by  exploiting  the  difference  between  the  instability  time  scale  and  the  period  of 
standing  waves.  The  perturbation  results  compare  very  favorably  with  the  exact  frequency-domain  theory 
of  Rott.  The  perturbation  results  are  further  simplified  by  introducing  a  short-stack  approximation  which 
is  nmumerically  much  simpler  and  only  slightly  less  accurate.  An  approximate  expression  fort  the  critical 
temperature  gradient  accounting  for  viscous  effects  and  other  design  features  is  also  derived.  In  addition  to 
the  fundamental  mode  of  a  tube  closed  at  both  ends,  the  theory  includes  higher  modes  as  well  as  open-end 
boundary  conditions. 

Karpov,  S.  &;  Prosperetti,  A.,  Nonlinear  Saturation  of  the  Thermoacoustic  Instability.  J. 

Acoust.  Soc.  Am.  107,  3130-3147,  2000 

ABSTRACT:  A  weakly  nonlinear  theory  of  the  thermoacoustic  instability  in  gas-filled  tubes  is  developed  in 
the  timer  domain  by  exploiting  the  difference  between  the  instability  time  scale  and  the  period  of  standing 
waves.  By  carrying  the  expansion  to  fourth  order  in  the  perturbation  parameter,  explicit  results  for  the 
initial  growth,  nonlinear  evolution,  and  final  saturation  are  obtained.  The  dependence  of  the  saturation 
amplitude  upon  the  temperature  difference  in  the  stack,  the  tube  geometry,  stack  plate  spacing,  Prandtl 
number,  and  other  parameters  is  illustrated. 

Karpov,  S.  &  Prosperetti,  A.,  A  Nonlinear  Model  of  Thermoacoustic  Devices.  J.  Acoust.  Soc. 
Am .,  accepted 

ABSTRACT:  This  paper  presents  a  nonlinear,  time-domain  model  of  thermoacoustic  devices  based  on  cross- 
sectional  averaged  equations.  Heat  transfer  perpendicular  to  the  device  axis  -  which  lies  at  the  core  of 
thermoacoustic  effects  -  is  modelled  in  a  novel  and  more  realistic  way.  Heat  conduction  in  the  solid  surfaces 
surrounding  the  fluid  medium  is  included.  Contrary  to  the  previous  versions  of  this  model  (Watanabe  et 
ah,  J.  Acoust.  Soc.  Am.  102,  3484,  1997),  the  present  version  does  not  require  artificial  damping  and  is 
numerically  robust.  The  model  performance  is  illustrated  on  several  examples:  a  prime  mover,  an  externally 
driven  thermoacoustic  refrigerator,  and  a  combined  prime  mover /refrigerator  system. 

Karpov,  S.,  Prosperetti,  A.,  fe  Ostrovsky,  L.,  Nonlinear  oscillations  of  bubble  layers.  J.  Acoust. 
Soc.  Am .,  accepted 

ABSTRACT:  Due  to  the  large  compressibility  of  gas  bubbles,  layers  of  a  bubbly  liquid  surrounded  by  pure 
liquid  exhibit  many  resonances  that  can  give  rise  to  a  strongly  non-linear  behavior  even  for  relatively  low- 
level  excitation.  In  an  earlier  paper  (Druzhinin  et  al.,  J.  Acoust.  Soc.  Am.  100,  3570,  1996)  it  was 
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pointed  out  that,  by  exciting  the  bubbly  layer  in  correspondence  of  two  resonant  modes,  so  chosen  that 
the  difference  frequency  also  corresponds  to  a  resonant  mode,  it  might  be  possible  to  achieve  an  efficient 
parametric  generation  of  a  low-frequency  signal.  The  earlier  work  made  use  of  a  simplified  model  for  the 
bubbly  liquid  that  ignored  the  dissipation  and  dispersion  introduced  by  the  bubbles.  Here  a  more  realistic 
description  of  the  bubble  behavior  is  used  to  study  the  nonlinear  oscillations  of  a  bubble  layer  under  both 
single-  and  dual-frequency  excitation.  It  is  found  that  a  difference-frequency  power  of  the  order  of  1%  can  be 
generated  with  incident  pressure  amplitudes  of  the  order  of  50  lcPa  or  so.  It  appears  that  similar  phenomena 
would  occur  in  other  systems,  such  as  porous  water-like  or  rubber-like  media. 

Karpov,  S.,  Nonlinear  Phenomena  in  Thermoacoustics  and  Bubbly  Liquids 

Doctoral  dissertation,  The  Johns  Hopkins  University,  May  2000. 


1.2  Personnel 

The  work  described  here  is  the  result  of  a  collaboartion  among  the  following  people: 

1.  A.  Prosperetti,  PI 

2.  L.  Ostrovsky,  Zel  Technology /NO  A  A  ETL,  Boulder,  CO 

3.  M.  Watanabe,  doctoral  student  (currently  associate  professor  at  Fukuoka  University,  Japan) 

4.  H.  Yuan,  post-doctoral  fellow 

5.  S.  Karpov,  doctoral  student  (currently  with  PTC,  Needham,  MA) 
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Chapter  2 


A  Simplified  Nonlinear  Model  of 
Thermoacoustic  Devices 

2.1  Introduction 

In  recent  years  there  has  been  a  renewal  of  interest  in  thermoacoustic  devices,  both  prime  movers  and  heat 
pumps  (for  reviews  see  The  foundations  of  the  linear  theory  were  firmly  established  in  a  well-known  series 
of  papers  by 

The  linear  analytical  theory  of  thermoacoustic  devices  has  reached  a  high  degree  of  maturity  thanks  to 
the  work  of  many  authors  (Rott  1969,  1976,  1980,  1983,  and  others;  Merkli  &  Thomann  1975;  Yazaki  et  al. 
1980,  1987;  Arnott  et  al.  1991,  1992,  1994,  1995;  Atchley  1992,  1994;  Atchley  &  Kuo  1994;  Atchley  et  al. 
1990a,  1990b,  1992;  Olson  &  Swift  1994;  Raspet  et  al.  1993;  Swift  1992;  Swift  &  Keolian  1993;  Wheatley  et 
al.  1983;  Karpov  &  Prosperetti  1998);  excellent  reviews  have  been  provided  by  Wheatley  (1986)  and  Swift 
(1988). 

Progress  in  the  direction  of  nonlinear  phenomena  has  been  slower  due  to  the  complexity  of  the  phenomena. 
Only  limited  analytical  results  are  available  (Gaitan  &  Atchley  1993;  Gopinath  et  al.  1998;  Bauwens  1996, 
1998;  Karpov  &  Prosperetti  2000) .  The  bulk  of  the  work  has  been  carried  out  with  the  aid  of  the  sophisticated 
numerical  techniques  necessary  to  tackle  this  multi-dimensional,  multi-scale  problem  (Cao  et  al.  1996; 
Worlikar  &  Knio  1996,  1998,  1999;  Worlikar  et  al.  1998,  Besnoin  &  Knio  2001). 

The  difficulty  of  the  nonlinear  problem  is  unfortunate,  because  this  is  the  situation  presenting  the  greatest 
practical  interest.  For  example,  the  very  fact  that  the  thermoacoustic  instability  in  a  prime  mover  ultimately 
saturates  at  a  finite  pressure  amplitude  is  an  obvious  consequence  of  nonlinearity.1  More  subtly,  as  pointed 
out  by  Worlikar  et  al.  (1998),  some  of  the  differences  between  linear  theory  and  experiment  encountered  by 
Atchley  et  al.  (1990a)  are  to  be  attributed  to  nonlinearity  already  at  drive  ratios  as  low  as  2%.  Furthermore, 
it  may  be  argued  that  the  effective  exploitation  of  thermoacoustic  systems  requires  operation  in  the  nonlinear 
regime.  The  strongly  nonlinear  resonators  described  by  Ilinskii  et  al.  (1998)  are  a  particularly  interesting 
possibility. 

We  have  therefore  thought  it  desirable  to  formulate  a  simplified  nonlinear  model  of  intermediate  com¬ 
plexity  between  the  linear,  frequency-domain  treatments  typical  of  the  linear  theory,  and  the  full  multi¬ 
dimensional  calculations  of  direct  numerical  simulations.  For  this  purpose,  we  have  used  a  method  very 
common  in  the  study  of  nonlinear  compressible  flows  and  shock  waves,  namely  the  integration  of  the  gov¬ 
erning  equations  over  the  cross  section  of  the  device  (see  e.g.  Crocco  1958;  Landau  and  Lifshitz  1959).  This 
procedure  has  led  us  to  a  quasi-one  dimensional  model  that,  although  approximate,  is  realistic  and  appears 
useful  to  further  an  understanding  of  thermoacoustic  devices  in  the  nonlinear  regime. 

1Heat  flow  limitations  in  the  heat  exchangers  can  also  affect,  the  saturation  level  quantitatively;  the  existence  of  the  saturation 
phenomenon  itself,  however,  is  a  nonlinear  effect. 
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The  averaging  of  the  equations  over  the  cross  section  of  the  device  has  the  effect  of  removing  a  detailed 
information  of  the  exchange  processes  of  momentum  and  energy  between  the  working  fluid  and  the  solid 
structure  in  contact  with  it.  These  processes  -  and  in  particular  the  second  one  -  lie  of  course  at  heart 
of  thermoacoustic  phenomena  and,  therefore,  it  is  essential  to  reintroduce  the  lost  information  into  the 
mathematical  model  with  a  sufficient  degree  of  realism.  This  is  a  crucial  point  discussed  in  detail  in  section 
2.4. 

As  a  result  of  our  work  over  the  past  years,  we  have  developed  a  general  time-domain,  fully  nonlinear 
model  of  thermoacoustic  devices  that  is  capable  of  predicting  the  response  of  such  systems  to  cross-sectional 
area  variations,  stack  position,  fluid  properties,  and  many  other  design  variables.  While  the  model  is  not 
exact  due  to  the  averaging  over  the  cross  section,  in  the  frequency  domain  and  for  small  amplitudes  it 
reduces  to  the  standard  linear  theory.  In  addition,  it  is  able  to  describe  nonlinear  phenomena  such  as  the 
finite-amplitude  saturation  of  the  thermoacoustic  instability  while,  at  the  same  time,  remaining  substantially 
simpler  than  a  truly  multidimensional  model.  We  thus  hope  that  our  model  can  be  useful  in  filling  the  current 
gap  in  the  levels  of  description  available  for  thermoacoustic  systems. 


2.2  Mathematical  Model 


Most  thermoacoustic  devices  consist  of  an  acoustic  resonator  containing  different  heat  transfer  components 
(stack,  heat  exchangers,  etc.)-  Typically  the  dimensions  along  the  direction  of  the  particle  displacement  - 
the  resonator  “axis”  -  is  much  longer  than  the  transverse  one  and  this  circumstance  suggests  the  basis  for 
our  approximation.  We  recast  the  governing  equations  in  an  integrated  form  over  the  cross-section  of  the 
device  thus  reducing  the  model  to  one  dimension  in  space  (along  the  tube  axis)  and  time.  Effects  taking 
place  in  the  orthogonal  directions  (friction,  heat  transfer,  etc.)  are  to  be  accounted  for  approximately  by  the 
introduction  of  suitable  terms  in  the  equations. 

Consider  a  thermoacoustic  device  consisting  of  a  duct  of  variable  area  S(x).  The  coordinate  x  is  taken 
along  the  axis  of  the  device  that  is  not  necessarily  straight.  Upon  integrating  the  equation  of  continuity  over 
the  volume  delimited  by  two  neighboring  cross  sections  S(x)  and  S(x  -f  dx)  we  find  the  well-known  relation 


d  <  p>  1  dS  <  pu> 
dt  +  S  dx 


where  p  is  the  gas  density,  u  the  velocity  in  the  x-direction,  and  the  angle  brackets  indicate  the  cross  sectional 
average, 

<...>=  -±-f  dS  (...).  (2.2) 

S{x)  JS{X) 

The  vanishing  of  the  exact  velocity  field  on  the  lateral  walls  of  the  duct  has  been  applied  to  obtain  (2.1). 
Similarly,  the  momentum  equation  in  the  x-direction  becomes 


d  <  pu>  1  9  /  f1  ^  v.  ->*  dS 

-^  +  sai{s<f,“>)  +  ~ar-  +  s{<p>-p}di 


V 


(f  •  n)a 


(2.3) 


Here  p  is  the  gas  pressure,  r  the  viscous  stress  tensor,  and  the  overline  denotes  the  average  over  the  “wetted 
perimeter”  £,  i.e.  the  fines  along  which  the  cross-section  S  is  cut  by  solid  boundaries, 


P  =  ^  fcPM 

where  V  is  the  length  of  C.  The  unit  normal  n  is  directed  into  the  fluid  region.  The  viscous  component  txx 
has  been  neglected  in  deriving  (2.3). 

Before  turning  to  the  energy  equation  we  introduce  an  assumption  widely  used  in  Gas  Dynamics  (see  e.g. 
Crocco  1958  for  a  discussion),  namely  that  the  fields  are  approximately  uniform  over  the  cross  section.2.  As 

2The  assumption  of  uniform  pressure  over  the  cross  section  has  already  been  used  in  Thermoacoustics  e.g.  by  Tominaga 
(1995) 
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a  consequence  of  the  cross  sectional  uniformity,  (p)  ~  p  and  we  may  disregard  correlation  terms  writing  the 
average  of  products  as  products  of  averages.  This  approximation  is  addressed  quantitatively  in  section  5. 
The  effects  of  non-uniformities,  such  as  wall  drag,  is  accounted  for  in  an  approximate  manner.  For  the  wall 
shear  stress  we  write 

— (r^n)x  =  V{pu) ;  (2.5) 

the  explicit  form  of  V(pu)  will  be  discussed  in  detail  in  the  next  section. 

With  these  assumptions,  and  the  omission  henceforth  of  explicit  indication  of  cross-sectional  averages, 
the  continuity  (2.1)  and  momentum  (2.3)  equations  become 


dp  1  d 

di  +  Sfa(SpU^ 


0, 


dpu  1  9  dp  x 

-dt+sdi{Spu')  +  te  -  ~v{pu)  • 

The  cross-sectional  averaging  procedure  applied  to  the  energy  equation  of  a  perfect  gas  gives 


d_ 

dt 


1  ,  1  2 

- 7P  +  o  PU 

7  —  1  2 


+ 


id_ 
S  dx 


uS 


7  ,  1  2 

- 7V  +  ~plt 

7  —  1  2 


=  5q  n> 


(2.6) 

(2.7) 


(2.8) 


where  q  is  the  heat  flux  and  7  is  the  ratio  of  specific  heats.  The  most  intense  heat  transfer  occurs  in 
the  directions  normal  to  x  and,  accordingly,  the  effect  of  axial  conduction  qx  has  been  disregarded.  This 
approximation  is  of  course  invalid  at  the  ends  of  the  tube.  We  return  on  this  point  in  subsection  2.1  below. 
Similar  to  (2.5),  we  introduce  two  heat  transfer  operators  H,  Q  by  writing 


V  BT 

-cpn  =  pcpU  (Tw  -  T)  -  -g^-Q(cppu) , 


(2.9) 


where  Tw(x,t )  is  the  (possibly  time-dependent)  temperature  of  the  solid  structure  (e.g.  the  stack  plates) 
at  x.  The  second  term  accounts  for  the  distortion  of  the  temperature  distribution  due  to  the  flow  in  the 
presence  of  a  mean  temperature  gradient.  This  ansatz  is  suggested  by  the  structure  of  the  temperature 
distribution  given  by  the  exact  linear  theory  and  will  be  discussed  in  detail  below  in  section  2.4. 

Now,  upon  using  (2.6)  and  (2.7)  to  eliminate  the  time  derivatives  of  p  and  u ,  the  energy  equation  (2.8) 
becomes3 

dT. 

CppH  ( Tw  -  T)  -  ~Q{pu)  +  uV(pu) 


dx 


(2.10) 


The  last  term  in  the  right-hand  side  represents  the  rate  of  conversion  of  mechanical  to  thermal  energy  by 
friction. 

The  mathematical  description  of  the  fluid  is  completed  by  the  specification  of  an  equation  of  state.  We 
assume  that  the  averaged  variables  are  related  by  the  perfect  gas  equation 


p  -  RpT , 


(2.11) 


where  R  is  the  universal  gas  constant  divided  by  the  gas  mass. 

Closure  of  the  model  requires  knowledge  of  the  temperature  distribution  Tw(x,  t)  along  the  axis  of  the 
device.  We  use  once  more  a  cross-sectional  average  over  the  solid  material  that  constitutes  the  stack. 
Although  it  is  not  necessary  to  specify  the  geometrical  structure  of  the  stack,  for  ease  of  exposition  we  will 
think  of  it  as  constituted  by  a  set  of  parallel  plates.  Area  averaging  is  justified  by  the  fact  that  the  plates 
are  usually  thin;  a  similar  remark  would  apply  to  pin  stacks  or  other  geometrical  arrangements.  In  this  way 
one  finds 


dTw 
psCs  dt 


BT 

pn(Tw-T)-^Q(pu) 


(2.12) 


3We  assume  that  cp  is  a  constant,  in  agreement  with  our  treatment  of  the  gas  as  perfect. 
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where  ps,  cs,  ks,  and  Ss  are  the  density,  specific  heat,  thermal  conductivity,  and  cross  sectional  area  of  the 
plates.  The  heat  exchangers  are  modelled  in  a  similar  way,  except  that  their  thermal  conductivity  is  taken 
to  be  infinite  so  that  their  temperature  is  spatially  uniform.4  For  the  sake  of  simplicity,  in  the  examples  that 
follow  we  assume  that  the  stack  plates  are  in  perfect  thermal  contact  with  the  cold  and  hot  heat  exchangers. 

Away  from  the  stack  region,  one  could  use  a  similar  formulation  for  the  resonator  tube  walls,  or  one 
could  couple  the  model  equations  with  the  conduction  equations  in  the  walls.  For  simplicity,  we  have  chosen 
to  assume  that  the  tube  wall  temperature  is  prescribed. 


2.3  Boundary  conditions 


A  variety  of  boundary  conditions  can  be  associated  with  the  mathematical  model  described  in  the  previous 
section  depending  on  the  situation  that  it  is  desired  to  model.  For  example,  for  the  prime  mover  case  with 
standing  waves,  it  might  be  reasonable  to  assume  that  the  end  walls  of  the  tube  are  rigid,  so  that  the  velocity 
vanishes: 

u  =  0  at  x  =  0,  x  —  L.  (2.13) 

The  momentum  equation  (2.7)  then  implies  that 

=  0  at  X  =  0,  x  =  L.  (2.14) 

OX 


Equations  (2.6)  and  (2.10)  written  at  the  endpoints  ( x  =  0,  L)  give  then 


dp  du 
dt  ^  dx 


=  0, 


dp 

dt 


Upon  eliminating  du/dx  one  finds 


du 


IP 


dx 


(7-l  )pcpH{Tw-T). 


dr 

dt 


7-1  T  dp 
7  p  dt 


+  H(TW-T). 


(2.15) 

(2.16) 


(2.17) 


This  relation  shows  that  a  knowledge  of  p  at  the  boundary  completely  specifies  T.  No  additional  boundary 
conditions  on  this  quantity  are  therefore  necessary.  Since  the  heat  transfer  term  H  is  very  small  outside 
the  stack  region,  this  relation  then  essentially  implies  the  adiabatic  pressure- temperature  relation  of  perfect 
gases. 

In  general,  the  relation  (2.17)  will  result  in  a  temperature  jump  from  the  T  given  by  this  relation  to 
the  end-wall  temperature.  This  is  a  consequence  of  the  neglect  of  axial  conduction  in  the  derivation  of  the 
energy  equation  (2.10).  Axial  conduction  would  introduce  a  term  d2T/dx 2,  important  only  near  the  end 
walls,  the  role  of  which  would  be  to  re-establish  continuity  of  temperature  by  means  of  a  thin  boundary 
layer.  The  temperature  in  this  layer  would  adjust  itself  so  as  to  match  the  value  given  by  (2.17).  This  is 
an  essentially  passive  process  with  negligible  effects  on  the  temperature  distribution  elsewhere  in  the  device 
and  can  therefore  be  disregarded.  If  desired,  however,  the  expression  for  H  can  be  adjusted  to  give  a  heat 
transfer  coefficient  at  the  tube  ends. 

For  the  case  of  standing  waves  forced  by  a  piston,  as  in  a  refrigerator  arrangement,  one  might  consider 
two  limit  cases.  The  piston  may  be  located  at  a  velocity  antinode,  corresponding  to  a  pressure  node.  In 
this  case  one  may  assume  a  vanishing  pressure  disturbance  and  an  imposed  fluid  velocity  at  this  location. 
Conversely,  for  a  piston  located  at  a  pressure  antinode,  one  may  use  an  imposed  pressure  condition  and  a 
zero  velocity  condition. 

4 It  appears  possible  to  modify  this  equation  to  account  for  possible  heat  flow  limitations  of  the  heat  exchangers  (e.g.,  by 
adding  a  term  mimicking,  by  means  of  a  heat  transfer  coefficient,  the  coupling  with  external  thermal  reservoirs),  although  have 
not  pursued  this  possibility. 
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These  models  can  of  course  be  refined.  For  example,  more  complicated  end-wall  impedances  may  be 
accounted  for  by  prescribing,  in  place  of  (2.13),  (2.14),  a  relation  coupling  pressure  and  velocity  at  x  =  0,  L. 
The  use  of  periodicity  conditions  at  the  tube  ends  would  permit  to  model  a  toroidal  travelling  wave  system, 
and  so  on. 


2.4  Exchange  Terms 


The  value  of  the  present  simplified  model  consists  in  its  ability  to  give  information  on  nonlinear,  time-domain 
processes  for  which  no  exact  theory  is  available  short  of  direct  multi-dimensional  numerical  simulations.  In 
order  to  accomplish  this  goal,  it  is  necessary  to  develop  a  suitable  time-domain  formulation  for  the  cross-axis 
exchange  terms  of  momentum,  V{u),  and  energy,  U(TW  -  T)  and  QdTwfdx .  For  this  purpose  we  turn  to 
the  exact  frequency- domain  theory  of  the  linear  regime. 

Swift  (1988,  Eq.  A4)  gives  the  following  expression  for  the  velocity  field  in  a  plane  channel  of  width  l: 


ui 


i  dpr 
cjpo  dx 


cosh  (1  +  i)y/Sy  \ 
cosh(l  +i)l/25v  ) 


where  y  is  measured  from  the  center  plane  of  the  channel  and  the  diffusion  length  8y  is  given  by 

Sv  = 

where  v  is  the  kinematic  viscosity.  From  this  equation  we  have 

du\ 


r  •  nx  =  -  p 


dy 


Vtanh(l+i) 


y=*/2 


UjpoSy  dx 


2  Sy 


The  mean  velocity  in  the  channel  (to  be  identified  with  our  u )  is 

re/2 

-if  2 

where 


=  -t 

tJ-i 


«i  dy  =  —  (1  ~  fv)-r~, 
6Jpo  dx 


fv  = 


tanh  (1  +  i)l/28y 
(l  +  i)£/28v 


From  these  two  expressions  we  find 


r(r  *  n)x  =  -iupo 


fv 


1  ~fv 


u . 


(2.18) 


(2.19) 


(2.20) 


(2.21) 


(2.22) 


(2.23) 


To  proceed  in  the  same  way  in  the  case  of  energy  exchange,  we  use  the  expression  for  the  temperature 
perturbation  in  the  channel  given  in  Eq.  (A10)  of  Swift  (1988): 


Tx  = 


cosh(l  +  i)y/SK 
pocp  [/  cosh  (1  4*  i)£/28k 


1  - 


1  dp1  dTu j 


a  cosh(l  +  i)y/8y  1  cosh  (1 i)y/ 8k 
a  —  l^cosh  (1  +  i)lj28y  a  —  1  cosh  (1  +  i)£/28K 


u)2po  dx  dx 

where  the  thermal  diffusion  length  8k  is  defined  similarly  to  8y  by 

(2a  1 


6k  = 


\f° 


Sy  , 


(2.24) 


(2.25) 
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in  which  a  is  the  thermal  diffusivity  and  a  the  Prandtl  number.  The  structure  of  this  result,  with  dependence 
on  both  Sv  and  5k,  shows  that  the  temperature  disturbance  is  determined  by  both  convection  and  conduction 
processes.  This  is  the  reason  for  our  ansatz  (2.9)  to  express  the  mean  energy  transport. 

Using  (2.24) ,  we  can  calculate  the  mean  of  Ti  over  the  channel  width,  to  be  identified  with  the  temperature 
disturbance  T  -  Tw  of  our  theory: 


T  —  Tw  =  (1  -fi<)—  + 


p0Cp  cj'2p0(l-a) 


[<t(1  -  fv)  -  (1  -  fn)\ 


dJ\, 

dx 


dpf 

dx 


(2.26) 


where,  as  in  (2.22), 


k  = 


tanh  (1  +  Vjlj  28k 


(1  +  i)(./28k 

Again  from  (2.24),  one  can  calculate  the  heat  flux  into  the  boundaries: 


_  U 

-q  n  =  7r 

Ots 


. ,  p'  fv  —  fie  9Tw 

ifi< - H 


PoCp  (1  -  fv)(  1  -  cr)  dx 

Upon  using  (2.21)  and  (2.26),  we  then  find 

V -  iuk  fK  .  k  dTw  (  1  a  \ 

-sq'“  =  — rrfc<T'T”)+S(T^)-8 


(2.27) 


(2.28) 


(2.29) 


The  functions  fv,  }k  defined  in  (2.22),  (2.27)  are  appropriate  for  a  plane  channel  geometry.  As  shown 
by  Rott  (1969),  for  a  circular  tube  with  radius  ro,  one  has 


,  _  ~  l)(ro /Sy)) 

Jv  Z(i-l)(ro/Jv)J0((*-l)(ro/^))  ’ 


(2.30) 


with  a  corresponding  expression  for  fu- 

The  expression  (2.28)  for  the  heat  transfer  term  may  appear  somewhat  unconventional  and  it  is  therefore 
appropriate  to  discuss  its  physical  significance.  It  is  easy  to  show  that  the  coefficient  of  Tw  —  T  in  the  first 
term  is  just  a  generalization  of  the  standard  heat  transfer  coefficient.  Indeed,  if  we  consider  the  limit  of 
steady  flow,  in  which  ui  — 0,  one  readily  finds 

qn  ~  8j(Tw-T),  (2.31) 

which  coincides  with  the  standard  result  for  steady  fully  developed  flow  in  a  channel  (see  e.g.  Incropera 
and  DeWitt  1996,  p.  430).  The  complex,  ^-dependent  form  that  appears  in  (2.28)  accounts  for  the  phase 
relation  that  is  crucial  in  an  oscillatory  flow.  The  same  parametrization  of  this  term  as  for  the  momentum 
transfer  coefficient  is  justified  by  the  similarity  between  the  two  transfer  processes,  which  is  also  apparent 
by  the  similar  functional  form  of  the  two  terms  (cf.  Eqns.  2.22  and  2.27). 

Of  greater  interest  is  the  second  term  in  (2.28)  which  does  not  arise  in  conventional  heat  transfer,  in 
which  wall  temperature  gradients  are  usually  not  accounted  for  explicitly.  The  physics  behind  this  term  may 
be  better  understood  if  it  is  noted  that  linearization  of  the  energy  equation  (2.10),  use  of  the  equation  of 
state  (2.11),  and  of  the  continuity  equation  (2.6)  permit  us  to  write 


pcP 


ituT*  +  (Td*  Q)  u 


dTw 

dx 


-  iwj/  -  -nr , 


(2.32) 


where  V  =  T  —  Tw,  p'  =  p  -  Po-  This  relation  shows  that  Q  represents  a  correction  to  the  convective 
transport  term  the  origin  of  which  lies  in  the  non-uniformity  of  the  velocity  distribution  over  the  channel. 
Consider  the  enthalpy  convected  by  the  fluid  during  a  time  dt  across  a  cross  section  located  at  x:  the  particles 
very  near  the  wall  will  carry  an  enthalpy  close  to  cpTw(x ),  because  they  will  have  moved  very  little  due  to 
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the  no-slip  condition.  However,  particles  further  way  from  the  wall  will  come  from  further  upstream,  where 
the  temperature  is  significantly  different  from  Tw{ x).  Thus,  the  enthalpy  convected  by  the  average  velocity 
u  should  not  be  cpTw(x),  but  should  be  modified:  the  factor  1  +  Q  has  precisely  the  role  to  effect  this 
modification.  This  analysis  is  supported  by  the  fact  that,  in  the  absence  of  viscosity  (a  =  0)  the  velocity- 
distribution  is  uniform  and  one  finds  indeed  Q  =  0  because,  then,  all  the  fluid  particles  transport  the  same 
enthalpy. 

Returning  now  to  the  modelling  of  the  exchange  terms  V,  U,  and  Q,  we  note  that,  upon  setting 

V{fm)  =  D{u)p0u,  %{TW  -  T)  =  H{u)  (Tw  -  T) ,  Q(pu)  =  Q(lj)  Pou ,  (2.33) 

our  formulation  will  match  the  exact  linear  theory  expressions  (2.23)  and  (2.29)  provided  that 


D{u)  =  iw 


fv_ 

1  —  fv 


(2.34) 


H{oS)  =  ico 


f_K 

W*  ’ 


CM 


1  /  fv  _  vfr<  \ 

1- a  V 1  -  fv  1  ~  fi<) 


(2.35) 


While  these  definitions  reproduce  the  linear  results,  they  are  unsuitable  for  the  present  purposes  of 
developing  a  time-domain  formulation.  In  our  earlier  work  (Yuan  et  al.  1997),  our  first  attempt  was  to 
simply  take 

V(pu)  =  D  ^1  +  pu ,  (2.36) 

with  D,Q&  constants  determined  in  such  a  way  that 


D  (1  +  =  jD(cji)  , 


(2.37) 


where  wi  is  the  real  part  of  the  frequency  of  the  fundamental  mode  of  the  system;  similar  forms  were 
postulated  for  the  other  exchange  terms.  This  choice  was  based  on  the  fact  that  the  fundamental  mode  at 
frequency  uq  is,  usually,  the  most  important  one  and  it  was  therefore  felt  desirable  to  simulate  it  as  precisely 
as  possible.  This  simple  choice  proved  however  inadequate  as  it  caused  many  high-order  modes  of  the  system 
to  become  unphysically  unstable,  which  caused  difficulties  in  numerical  work. 

We  have  therefore  decided  to  improve  on  this  approach  by  following  a  suggestion  of  Achard  &  Lespinard 
(1981),  who  were  interested  in  finding  a  quasi-one  dimensional  formulation  for  the  time-dependent  viscous 
flow  of  a  fluid  in  a  duct.  While,  with  the  assumption  of  fully  developed  flow,  the  problem  is  linear  and 
can  be  solved  by  Laplace  transform  methods,  they  realized  that  no  form  for  the  drag  term  short  of  a  time 
convolution  can  capture  exactly  the  physics  of  the  process.  For  this  reason,  they  proposed  an  approximation 
in  which  the  momentum  transfer  term  was  not  given  by  a  single  constant,  but  was  found  by  solving  an 
ordinary  differential  equation  in  time.  The  simplest  possibility  beyond  (2.36)  is  to  write 

aD^  +  V  =  bD  (l  +  cD^j  pu ,  (2.38) 

with  aD,bo,CD  three  suitable  real  constants.  If  ao  =  0,  this  relation  reduces  to  the  earlier  choice  (2.36) 
but,  in  general,  it  contains  one  additional  parameter  that  can  be  chosen  so  as  to  improve  the  representation 
of  the  fluid-solid  momentum  transfer.  We  proceed  in  a  similar  fashion  for  the  other  exchange  terms  as  well: 

aH W  +  n  =  bH  (x  +  CH§i)  {Tw  ~  T)  ’  (2'39) 
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ReQ,  Im Q 


In  the  frequency  domain  all  these  prescriptions  reduce  to  (2.33)  with 


D,  H,  Q  replaced  by 


D{(jj)  —  bo 


1  4*  iuJCjy 
1  +  iU)CL£> 


(2.41) 


etc.  We  are  now  at  liberty  to  impose  three  conditions  in  order  to  determine  the  model  parameters  a,  b,  c.  As 
before,  in  view  of  the  importance  of  the  fundamental  mode  at  frequency  wi,  we  impose  that 


D(u>i)  =  bo 


1  +  iuicp 

1  +  iujian 


D(uj i) , 


(2.42) 


with  D(ui)  given  by  (2.34);  here,  as  before,  is  the  real  part  of  the  complex  eigenfrequency  uq  of  the  first 
mode.  Upon  separating  real  and  imaginary  parts,  this  relation  gives  two  equations. 

For  the  remaining  condition  there  is  a  considerable  latitude.  Physically,  since  the  transverse  dimen¬ 
sions  of  the  flow  passages  in  the  stack  are  usually  much  smaller  than  the  stack  length,  one  can  make  the 
approximation  of  fully  developed  flow  in  the  stack  which  has  the  effect  of  rendering  the  left-hand  side  of 
the  exact  conservation  equations  linear.  The  nonlinear  mode-mode  coupling  is  essentially  confined  to  the 
region  outside  the  stack,  and  is  therefore  described  by  the  left-hand  sides  of  the  exact  conservation  equations, 
which  are  well  approximated  by  the  left-hand  sides  of  (2.6)  to  (2.10).  The  process  that  the  exchange  terms 
must  model,  therefore,  is  essentially  a  linear  process,  which  suggests  to  try  to  proceed  on  a  mode-by-mode 
basis.  Ideally,  therefore,  the  best  approach  would  be  to  choose  the  remaining  constant  so  as  to  minimize  the 
difference  between  the  exact  and  approximate  linear  spectra.  This  objective  can  be  achieved,  but  once  more 
at  the  cost  of  solving  the  eigenvalue  problem  exactly. 

Hence,  we  propose  a  simpler  alternative.  In  a  thermoacoustic  prime  mover,  it  is  usually  the  first  mode  that 
is  unstable.  As  its  amplitude  grows,  it  loses  energy  to  the  second  and  higher  modes  by  non-linear  couplings. 
The  second  mode  is  excited  the  most  and,  since  it  is  stable,  it  will  represent  the  greatest  energy  sink  for 
the  system.  The  situation  is  similar  in  a  thermoacoustic  refrigerator  where  the  forcing  typically  energizes 
the  fundamental  mode  the  most  with,  again,  the  second  mode  providing  the  greatest  energy  loss  (after  the 
fundamental).  These  considerations  suggest  that  an  effective  second  condition  for  the  determination  of  the 
last  free  constant  in  the  exchange  terms  is  to  impose  that  the  damping  of  the  second  mode  be  correctly 
described.  We  cannot  use  condition  (2.42)  for  the  second  mode  since  the  real  and  imaginary  parts  of  the 
relation  (2.42)  give  two  equations,  but  only  one  free  constant  remains  available.  To  choose  the  right  condition 
we  use  as  a  guide  an  earlier  result  (Karpov  and  Prosperetti  1998  Eqs.  33,  34;  Karpov  and  Prosperetti  2000  Eq. 
5.4)  according  to  which  the  linear  growth  (or  damping)  rate  of  the  generic  n-th  mode  can  be  approximated 

by 


Im  <2>n 


2Vpfar 


f 


dxS 


^  ReD(u}n) 
Un 


( dPn\ 2 

V  dx  ) 


+  7wn  Ke//(can)  Pn  (Pn  7ZRnTw)  4-  (7  1)  ImQ(uVi)  Cp 


dT^p  dPn 

dx  ”  dx 


(2.43) 


This  result  shows  that  Im  Qin  depends  on  Re  Re  H(ton),  and  Im  <2(u>„),  which  suggests  the  following 

conditions: 

ReD(uj2)  =  bD  Ren  +tU)2C^,  (2.44) 

1  +  iw-2  an 

R eif(w2)  =  bH  Re  1  +  ,  Im<3(w2)  =  bQ  Im  1  *  luj2C9  .  (2.45) 

V  '  1  +  tU2aH  1+IU2CLQ 

Here  us?  is  the  real  part  of  the  complex  eigenfrequency  of  the  second  mode  In  order  to  demonstrate  the 
quality  of  the  approximations  thus  obtained,  in  Figs.  2.1  we  compare  the  ^-dependence  of  the  functions  D 

5It  is  seen  from  (2.35)  that  Q  is  a  combination  of  D  and  H.  However,  the  presence  of  the  factor  iu)  in  these  two  quantities 
shows  that,  the  physical  effect  accounted  for  by  D  and  H  is  in  phase  with  the  particle  displacement,  rather  than  the  velocity, 
as  the  effect  of  Q.  For  this  reason,  in  the  time  domain,  it  is  not  desirable  to  derive  the  expression  of  Q  from  those  of  V  and  71, 
but  to  proceed  directly  as  in  (2.35). 
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X  A  =  0.909 

o' 


Figure  2.2:  Damping  constant  of  the  first  20  eigenmodes  for  the  system  and  conditions  described  in  the  text; 
the  stack  is  located  at  xs/L  =  0.909.  The  circles  are  the  exact  results  of  the  linear  theory.  The  squares  and 
triangles  correspond  to  the  differential  equation  approximation  (2.38)  to  (2.40)  to  the  momentum  and  energy 
exchange  terms  and  to  the  simple  explicit  form  (2.37)  respectively.  The  lines  are  only  meant  as  guides  to 
the  eye. 


and  Q  given  by  (2.38),  (2.40)  with  the  exact  expressions  D,  Q  given  in  (2.34),  (2.35).  The  Prandtl  number 
is  cr  =  0.71  and,  at  the  frequency  oq,  <5A-/f  =  0.34.  Since  the  functional  dependence  of  V.  on  u  is  the  same 
as  that  of  V,  we  do  not  show  results  for  this  quantity. 

In  a  nonlinear  setting,  the  partial  time  derivative  in  (2.38),  (2.39),  and  (2.40)  could  be  replaced  by  the 
convective  derivative  d/dt  +  ud/dx.  We  have  taken  this  approach  in  Yuan  et  al.  (1997),  but  found  that 
the  difference  is  minimal,  probably  because  the  flow  in  the  most  of  the  stack  is  essentially  fully  developed. 
Hence,  we  retain  the  form  (2.38),  (2.39),  and  (2.40)  in  the  nonlinear  version  of  the  model  as  well. 

As  formulated,  it  would  seem  that  the  determination  of  the  model  requires  the  solution  of  the  exact  linear 
eigenvalue  problem,  at  least  for  the  determination  of  uq  and  Actually  this  task  can  be  considerably 
simplified  on  the  basis  of  the  following  observation.  In  principle,  the  linear  eigenfrequencies  of  the  system 
are  to  be  found  by  solving  the  linearized  version  of  the  model,  which  can  be  shown  to  reduce  to  (Watanabe 
et  al.  1997): 


ld_ 
S  dx 


Sc\  dp 
1  +D/(iu)  dx 


H 

+  uj-p  +  — 
IU 


+ 


C 7  l)cp  dd' w  ^dp 
1+D/M  dx  ^  dx 


S  dp 
1  +  D/(iuj)dx 


+  l^P  i 


(2,46) 


where  c2A  =  ~f7lTw  is  the  local  adiabatic  sound  speed  and  proportionality  to  exp  (iujt)  has  been  assumed.  It 
is  found  numerically  that  the  real  part  of  the  eigenfrequencies  is  very  well  approximated  by  neglecting  the 
exchange  terms  in  the  eigenvalue  equation  (2.46)  to  find 


(2.47) 


1-15 


A  further  simplification  is  obtained  by  using,  for  ca ,  a  value  based  on  the  average  temperature  of  the 
resonator: 

rL 

(2.48) 


Te  —  ^  J  Tw  dx. 

In  particular,  for  a  cylindrical  tube  closed  at  the  two  ends,  one  may  take 


a?i  ~  —  ^//y7ZTe  ,  lo 2  ~  2cji  . 


(2.49) 


Approximations  are  also  available  for  tubes  of  non-uniform  cross  section  (see  e.g.  Rayleigh  1896),  although 

(2.49)  may  give  an  adequate  approximation  in  many  of  these  cases  as  well. 

As  examples  of  the  effect  of  the  previous  approximations  on  the  linear  spectrum  of  the  problem,  Figs.  2.2, 
2.3,  and  2.4  show  the  damping  constant  of  the  first  20  eigenmodes  as  a  function  of  the  mode  number  for 
three  different  positions  of  the  stack  midpoint,  xs/L~ 0.91,  0.77,  0.68;  instability  corresponds  to  a  negative 
value  of  the  quantity  shown.  The  circles  connected  by  the  solid  line  are  the  exact  result  of  the  standard 
linear  theory  (2.46).  Note  that  in  all  the  examples  only  the  first  mode  is  unstable.  The  results  shown  by 
the  squares  and  the  dashed  line  correspond  to  the  differential  equation  formulation  (2.38)  to  (2.40),  while 
the  triangles  connected  by  the  dotted  line  correspond  to  the  simple  formulation  (2.34)  and  (2.35)  for  the 
exchange  operators.  Although  both  models  give  good  results  for  the  first  few  modes,  the  simple  formulation 
(2.34)  to  (2.35)  makes  some  higher  modes  negatively  damped,  i.e.,  unstable;  this  circumstance  was  at  the 
root  of  the  numerical  difficulties  encountered  in  our  previous  work  (Yuan  et  al.  1997).  On  the  other  hand, 
the  differential  equation  approximation  (2.38)  to  (2.40)  gives  a  good  agreement  with  the  exact  results  for 
the  first  few  modes  and  reasonable  estimates  for  the  higher  modes.  Since  the  fraction  of  energy  contained 
in  these  higher  modes  is  usually  very  small,  the  associated  error  is  likely  to  be  acceptable. 

These  numerical  results  are  for  a  system  consisting  of  a  1  m-long  cylindrical  tube  filled  with  helium 
(7=5/3,  cr—0.71,  cp= 5.2  kJ/kg  K).  The  stack  has  a  length  of  30  mm  and  consists  of  parallel  plates  of 
negligible  thickness  spaced  by  0.77  mm.  The  radius  of  the  tube  is  assumed  to  be  large  so  that  drag  and 
heat  transfer  effects  are  small  outside  the  stack  region.  The  mean  pressure  is  307  kPa.  The  tube  wall 
temperature  at  the  left  of  the  stack  is  Tc  =  293  K,  and  at  the  right  Tu  =  415  K  so  that  T#  —  Tq  —  122  K. 
In  the  stack  region  the  wall  temperature  is  independent  of  time  and  a  linear  function  of  position.  Over  the 
temperature  range  of  interest  the  thermal  conductivity  data  was  fitted  by  a  linear  function  of  temperature 
as  k  —  0.151  +  3.228  x  10“4(T  -  300),  with  k  in  W/m  K  and  T  in  K.  We  have  included  this  effect  in  the 
calculations  as  the  value  of  k  determines  the  thermal  penetration  thickness,  and  has  therefore  a  significant 
impact  on  the  heat  transfer  parameters.  The  frequencies  uq  and  U2  were  evaluated  using  the  approximation 

(2.49) . 

The  present  model  can  accommodate  tubes  with  a  non-uniform  cross  section.  To  illustrate  this  effect,  we 
show  in  Fig.  2.5  results  similar  to  Fig.  2.4  but  including  the  blockage  caused  by  a  finite  stack  plate  thickness 
with  S stack /Stube  =  0.75.  As  another  example,  in  Figs.  2.6  and  2.7  we  consider  a  tube  with  a  cross-sectional 
area  given  by 

1 


S(x) 

5(0) 


=  < 


[l  +  Z  cos2{5tt/2  (2 x/L  -  l)}f 


0<|  L 
$L  <  x  <  g L 
|  L  <x  <L 


(2.50) 


for  Z  equal  to  0.4  and  -0.4  respectively;  the  stack  position  is  xs/L= 0.684  and  the  thickness  of  the  stack  plates 
is  neglected.  In  all  these  examples  the  frequencies  c*q  and  uq  were  evaluated  using  the  simple  approximation 

(2.49) ,  which  remains  fairly  accurate  in  spite  of  the  cross  section  variation.  For  example,  for  the  case  of 

(2.50)  with  Z  =  0.4,  the  error  for  the  fundamental  mode  is  less  than  3%  and  for  the  second  mode  about 

10%. 
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xs/L  =  0.774 


Figure  2.3:  Damping  constant  of  the  first  20  eigenmodes  for  the  system  and  conditions  described  in  the  text; 
the  stack  is  located  at  xs/L  =  0.774.  The  circles  are  the  exact  results  of  the  linear  theory.  The  squares  and 
triangles  correspond  to  the  differential  equation  approximation  (2.38)  to  (2.40)  to  the  momentum  and  energy 
exchange  terms  and  to  the  simple  explicit  form  (2.37)  respectively.  The  lines  are  only  meant  as  guides  to 


xs/L  =  0.684 


Figure  2.4:  Damping  constant  of  the  first  20  eigenmodes  for  the  system  and  conditions  described  in  the  text; 
the  stack  is  located  at  xs/L  =  0.684.  The  circles  are  the  exact  results  of  the  linear  theory.  The  squares  and 
triangles  correspond  to  the  differential  equation  approximation  (2.38)  to  (2.40)  to  the  momentum  and  energy 
exchange  terms  and  to  the  simple  explicit  form  (2.37j)7respectively.  The  lines  are  only  meant  as  guides  to 
the  e3^e. 


stack/ tube 


Figure  2.5:  Damping  constant  of  the  first  20  eigenmodes  for  the  same  case  as  in  the  previous  figure,  but 
including  the  25%  blockage  (75%  porosity)  of  finite-thickness  plates. 


2.5  Numerical  Method 

With  the  expressions  (2.36),  (2.39),  and  (2.40),  and  the  definitions 

m  =  pu , 


e  = 


1  .  1  2 

- 7'P  +  ~PU  . 

7  —  1  2 


(2.51) 


for  the  momentum  flux  m  and  total  energy  e,  it  is  easy  to  verify  that  the  continuity,  momentum,  and  energy 
equations  (2.6),  (2.7),  and  (2.10)  may  equivalently  be  written  as 


dp  dm  m  dS  _  ^ 
dt  dx  S  dx  ’ 


dm  d  .  .  mu  dS  n  N 

+  _(mu  +  p)  +  __  =  0  =  -V(m), 


dt  dx 


de  +2-  [U(e  +/>)]  +  ^  (e  +  p)g  =  pcpH(Tw  -  T)  -  cP^Q(m) . 


dt  "  dx  l  K  '  5 

We  rewrite  this  system  of  equations  in  vector  form  as 

<9w  dF 


(2.52) 

(2.53) 

(2.54) 


dt+dx+s  b’ 


(2.55) 


where  w  and  the  flux  vector  F  are  given  by 


* 

p 

m 

w  = 

m 

,  F  = 

mu  4-  p 

* 

e 

(e  +  p)  u 

(2.56) 
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Figure  2.6:  Damping  constant  of  the  first  20  eigenmodes  for  a  resonator  with  a  wider  midsection  according 
to  (2.50)  with  Z  =  0.4;  see  text  for  a  description  of  the  system  modelled.  The  stack  is  located  at  xs/L  = 
0.684  as  in  the  previous  two  figures.  The  circles  are  the  exact  results  of  the  linear  theory.  The  squares  and 
triangles  correspond  to  the  differential  equation  approximation  (2.38)  to  (2.40)  to  the  momentum  and  energy 
exchange  terms  and  to  the  simple  explicit  form  (2.37)  respectively.  The  lines  are  only  meant  as  guides  to 
the  e}'e. 


Figure  2.7:  As  in  the  preceding  figure,  but  with  Z  —  -0.4. 
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while  the  vector  s  accounting  for  the  effect  of  changes  in  the  cross-sectional  area  and  the  forcing  b  are  given 

by  .  , 


m 

0 

s  = 

mu 

i  ds 

s  dx  ’ 

b  = 

—V(m) 

( e+p)u 

pcpH(Tw  -  T)  -  cp(dTw /dx)Q(m) 

Our  first  attempt  at  solving  Eq.  (2.55)  was  based  on  a  straightforward  centered-difference  spatial  dis¬ 
cretization  with  a  predictor-corrector  time  stepping  procedure  (Prosperetti  &  Watanabe  1994).  We  found 
that,  whenever  conditions  were  such  that  quasi-shocks  developed,  a  series  of  spurious  grid-dependent  oscil¬ 
lations  also  appeared.  Such  oscillations  are  a  well-known  numerical  artifact  affecting  computations  in  the 
presence  of  steep  gradients  (see  e.g.  Roe  1986;  Fletcher  1988),  and  their  elimination  has  motivated  a  large 
amount  of  research.  While  a  complete  bibliography  would  be  out  of  place  here,  it  may  be  useful  to  cite  the 
review  by  Roe  (1986)  and  a  few  other  papers  (Sod  1978;  Harten  1983;  Osher  1984;  Osher  and  Chakravarthy 
1984;  Sweby  1984;  Harten  et  al.  1986;  Harten  and  Osher  1987).  This  effort  has  led  to  a  new  family  of 
schemes  for  hyperbolic  equations  known  as  Total  Variation  Diminishing  (TVD)  schemes.  The  name  is  a 
consequence  of  the  definition  of  the  total  variation  TV(un)  of  a  grid  function  {uf },  i  =  1, 2, . . . ,  N  +  1  at 
time  tn: 

TV(uB)  =  f>r+1-<l,  (2-58) 

i~  1 

and  of  the  fact  that  these  schemes  have  the  property  that  TV(un)  is  a  non-increasing  function  of  time: 

TV(un)  -  TV(un+1).  (2.59) 

Evidently,  a  TVD  scheme  cannot  produce  an  oscillatory  solution  starting  from  monotonic  initial  data.  We 
have  found  that  the  scheme  proposed  by  Harten  (1983)  proved  suitable  for  our  problem. 

The  system  (2.55)  is  first  discretized  explicitly  in  time  as 


w 


"+1-w"  F»1/2-P?1/2 


At 


~  + 


i+ 1/2 i- 1/2 

Ax 


+  S 


(2.60) 


where  superscripts  indicate  time  levels  and  subscripts  spatial  nodes.  The  essential  aspect  of  the  method  is 
the  manner  in  which  the  modified  fluxes  F  are  specified  in  terms  of  the  auxiliary  quantities 


(2.61) 
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(2.62) 

(2.63) 

(2.64) 

(2.65) 
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where  the  quantities  carrying  a  half-integer  subscript  are  calculated  as  arithmetic  averages,  e.g.  u  —  ki^i+i  4- 
U{).  Then 


=  \  (F«  +  F.)  +  ^  E  [s!°  +  -  <3(-!?1/a  +  if )  a'0]  R';,/a  ,  (2-66) 


where 


M)  = 


At 
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(«)  _  gj+i  9i 
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a ■ 


(2.67) 


The  quantities  g\e)  appearing  in  these  definitions  are  a  correction  to  the  components  of  the  flux  F  along  the 
characteristic  directions  and  are  the  smaller  one  in  modulus  between 


[<9(^+i/2)  -  ^1/2]  “I+1/2  -  and  [^(^-1/2)  -  ^-1/2]  “1-1/2  • 


(2.68) 


This  flux  correction  is  introduced  to  account  for  the  discretization  error  and  guarantees  second-order  accuracy 
in  space. 

The  last  quantity  to  be  defined  is  the  function  Q(x)  that  may  be  considered  as  a  modified  |x|.  Specifically, 
following  Harten  (1983),  we  take 


Q( 


*>  "  {  12 


2/(4e)-f-c  for  \x\  <  2e 


for 


Id  >  2e 


(2.69) 


with  e  =  0.1.  This  quantity  plays  the  role  of  an  artificial  viscosity. 

The  other  terms  in  (2.60)  are  treated  in  a  more  conventional  way.  The  components  of  b  are  written  as 
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(2.71) 


The  spatial  derivatives  appearing  s  are  discretized  by  central  differences. 

After  the  quantities  p,m,  and  e  at  time  step  (n  +  1)  are  calculated,  the  gas  temperature  is  obtained 
from  the  equation  of  state  (2.11)  and  the  wall  temperature  T”+1  is  found  from  a  discretization  of  the  heat 
conduction  equation  (2.12): 
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The  last  step  is  to  update  the  exchange  terms: 
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t~l  - 
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The  values  of  viscosity  and  thermal  diffusivity  are  evaluated  at  the  local  instantaneous  prevailing  conditions. 

Since  the  integration  is  explicit  in  time,  the  preceding  formulae  are  sufficient  to  construct  the  solution 
at  all  the  interior  nodes  at  time  level  tn+1  starting  from  the  known  values  at  tn.  The  solution  at  the  two 
boundary  nodes  is  calculated  from  the  boundary  conditions  appropriate  for  each  particular  problem. 

The  time  step  At  is  adjusted  so  that  the  maximum  Courant  number  At|a(i)|/Ax  (with  the  all))s  defined 
in  2.62)  remains  below  0.4. 

It  is  well  known  (see  e.g.  Cao  et  al.  1996;  Worlikar  &  Knio  1999)  that  sharp  temperature  gradients  exist 
near  the  ends  of  the  stack,  while  the  fields  are  smooth  away  from  the  stack  region.  For  this  reason,  the  use 
of  a  variable  spatial  node  spacing  improves  efficiency. 


2.6  Results 

In  order  to  illustrate  the  behavior  of  the  model  described  before  we  now  consider  several  examples:  a  prime 
mover  with  a  temperature  gradient  above  onset,  an  externally  driven  thermoacoustic  refrigerator,  and  a 
prime  mover /refrigerator  combination. 

2.6.1  A  Prime  Mover 

The  first  example  is  a  model  of  the  prime  mover  system  studied  by  Atchley  et  al.  (1990b)  and  Atchley  (1992, 
1994).  It  consists  of  a  38.2  mm-diameter  tube  with  a  length  of  99.89  cm,  a  stack  of  35  stainless  steel  plates 
located  90.13  cm  from  the  cold  end,  and  two  heat  exchangers.  The  plates  are  35  mm-long,  with  a  thickness 
of  0.28  mm  and  a  spacing  of  0.77  mm.  We  take  handbook  values  for  the  physical  properties:  ps  =  7,900 
kg/m3,  cs  =  480  J/(kg  K),  ks  =  14.9  W/(m  K).  The  cold  heat  exchanger  at  the  left  end  of  the  stack  consists 
of  two  identical  structures  separated  by  1.5  mm,  each  with  25  nickel  plates  0.45  mm  thick,  spaced  by  1.04 
mm  and  10.2  mm  long.  The  hot  heat  exchanger  is  attached  to  the  right  of  the  stack  and  is  built  as  the 
cold  one  except  that  it  consists  of  only  one  7.62  mm-long  structure.  The  area  blockage  is  about  30%  (70% 
porosity)  in  the  heat  exchangers  and  27%  (73%  porosity)  in  the  stack. 

A  feature  to  be  noted  in  this  system  is  the  difference  between  the  number  of  stack  and  heat  exchanger 
plates.  For  the  conditions  of  the  experiment  the  thermal  penetration  length  at  the  temperature  of  the  cold 
heat  exchanger  is  approximately  0.19  mm,  which  is  about  25%  of  the  stack  plate  spacing  but  only  18%  of  the 
heat  exchanger  plate  spacing.  At  the  hot  heat  exchanger  temperature,  5k  ~  0.38  mm  and  the  corresponding 
ratios  are  49%  and  37%.  Since  a  gas  particle  only  exchanges  heat  if  it  is  within  1-2  thermal  penetration 
lengths  from  the  plates,  as  sketched  in  Fig.  2.8  it  may  expected  that  some  of  the  stack  plates  are  effectively 
not  contributing  as  the  gas  particles  that  move  along  them  make  only  an  imperfect  thermal  contact  with 
the  heat  exchangers,  particularly  at  the  cold  end  where  6k  is  smaller. 

To  account  in  a  rough  way  for  this  effect  one  may  say  that  only  as  many  stack  plates  as  there  are  heat 
exchanger  plates  take  part  in  the  heat  transfer.  Alternatively,  the  effect  of  each  stack  flow  passage,  as 
represented  by  the  heat  exchange  terms  K  and  Q,  should  be  reduced  by  a  factor  Kh  =  25/35  ~  0.71. 

In  order  to  show  that  this  is  a  reasonable  estimate  we  first  consider  the  linear  case.  With  the  correction 
factor  I<h ,  upon  linearization  and  in  the  frequency  domain,  in  place  of  (2.46)  the  system  (2.6)  to  (2.10)  leads 
to  the  following  eigenvalue  equation  for  the  pressure  perturbation  p': 
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Heat  exchanger  plates 


Stack  plates 


Figure  2.8:  Schematic  representation  of  a  case  in  which  the  gap  between  stack  plates  is  half  that  between  heat 
exchanger  plates.  The  gas  particles  moving  between  points  A  and  B  exchange  heat  with  the  heat  exchanger 
plates  since  the  point  A  is  within  a  thermal  penetration  depth  5x<-  In  contrast,  the  gas  particles  moving 
between  points  C  and  D  do  not  exchange  heat  with  the  heat  exchanger  plates.  Therefore,  the  middle  stack 
plate  does  not  participate  in  the  heat  transport  along  the  stack. 
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(2.76) 


where,  as  before,  ca  is  the  local  adiabatic  sound  speed  and  D(co),  H(cj ),  and  Q(uj)  are  defined  in  (2.34)  and 
(2.35). 

In  the  experiments  of  Atchley  et  al.  (1990b)  the  temperature  difference  of  325  K  was  slightly  above 
onset.  A  solution  of  equation  (2.76)  gives  an  onset  temperature  difference  of  279.6  K  for  Kjt  =  1  and  319.1 
K  for  Kh  =  0.71.  From  the  data  reported  in  Atchley  1994,  difference  along  the  stack  of  379  K,  the  linear 
temporal  growth  rate  of  the  perturbation  is  5.0  s-1.  The  solution  of  Eq.  (2.76)  gives  a  linear  growth  rate  of 
11.87  s-1  for  Kh  —  1  and  5.64  for  Kh  =  0.71.  Obviously  the  correction  embodied  in  Kh  gives  a  much  better 
agreement  of  the  linear  theory  with  the  experimental  results  and,  for  this  reason,  we  feel  justified  in  using 
the  same  correction  factor  Kh  for  the  nonlinear  problem  as  well. 

Figures  2.9  and  2.10  show  the  transient  and  steady-state  temporal  waveforms  of  the  pressure  at  the 
cold  end  of  the  tube  for  the  case  with  a  temperature  difference  Th  -  Tc  =  368  K  described  by  Atchley 
et  al.  (1990b).  The  calculation  is  started  with  a  linear  temperature  distribution  in  the  £tack  and  a  small 
amplitude  of  the  first  system  normal  mode.  The  helium  mean  pressure  is  307  kPa.  Figure  2.11  shows  the 
temperature  deviation  from  the  initial  value  as  a  function  of  time  at  the  center  of  the  stack  and  1  mm  away 
from  its  cold  and  hot  ends.  The  temperature  near  the  cold  end  increases  and  that  near  the  hot  end  decreases 
rapidly  during  the  first  second,  after  which  they  slowly  reach  the  steady  state  values.  The  difference  between 
the  initial  and  steady  state  temperatures  is  approximately  5  K  and  6  K  for  the  cold  and  hot  ends  of  the 
stack,  respectively.  In  contrast,  the  temperature  in  the  middle  of  the  stack  decreases  slowly  to  its  steady 
state  value  (not  yet  achieved  in  this  figure).  The  final  stack  temperature  distribution  averaged  over  one 
cycle  is  shown  in  Fig.  2.12.  The  solid  line  is  the  mean  wall  temperature  and  the  dotted  line  the  mean 
gas  temperature.  One  can  see  that  there  are  near-jumps  in  the  temperature  of  the  solid  structure  at  the 
ends  of  the  stack.  The  correctness  (at  least  qualitative)  of  this  result  is  confirmed  by  the  experiments  and 
analysis  of  Brewster  et  al.  1997  and  the  calculations  of  Worlikar  (1997)  and  Worlikar  et  al.  (1998).  These 
near-jumps  effectively  reduce  the  available  temperature  difference  in  the  stack  thus  reducing  the  steady  state 
oscillation  amplitude  in  comparison  with  what  would  be  obtained  by  assuming  a  time-independent  linear 
stack  temperature  distribution  equal  to  Th  —  Tc-  This  effect  is  shown  in  Fig.  2.10  by  the  dotted  line  which 
represents  the  results  with  the  assumption  that  the  wall  temperature  in  the  stack  region  is  a  time-independent 
linear  function  of  position.  For  this  particular  case  the  difference  between  the  two  results  is  however  small. 
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Figure  2.9:  Transient  behavior  of  the  pressure  at  the  cold  end  of  the  tube  for  a  prime  mover  with  a  temper¬ 
ature  difference  Th  —  Tc  =  368  K  described  by  Atchley  et  al.  (1990b). 
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Figure  2.10:  Steady-state  temporal  waveform  of  the  pressure  at  the  cold  end  of  the  tube  for  the  prime  mover 
of  the  previous  figure.The  dotted  line  corresponds  to  a  fixed  linear  temperature  distribution  in  the  stack;  the 
solid  line  is  calculated  with  the  numerically  determined  steady-state  temperature  distribution  at  the  end  of 
the  transient. 
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Figure  2.11:  Temperature  deviation  from  the  initial  value  as  a  function  of  time  at  the  center  of  the  stack 
and  1  mm  away  from  the  cold  and  hot  ends  for  the  prime  mover  of  the  previous  two  figures.  Initially  the 
temperature  near  the  cold  end  increases  and  that  near  the  hot  end  decreases  very  rapidly.  The  difference 
between  the  initial  and  steady  state  temperature  is  approximately  5  K  and  6  K  for  the  cold  and  hot  ends  of 
the  stack,  respectively.  In  contrast,  the  temperature  in  the  middle  of  the  stack  decreases  slowly  to  its  steady 
state  value  (not  shown  in  this  figure).  The  final  stack  temperature  distribution  is  shown  in  the  next  figure. 


Figure  2.12:  The  final  mean  temperature  distribution  in  the  stack  (solid  line)  and  in  the  gas  (dotted  fine) 
for  the  prime  mover  of  the  previous  three  figures.  Note  the  jumps  in  the  solid  structure  temperature  at  the 
ends  of  the  stack.  T25 


Figure  2.10  should  be  compared  with  Fig.  4  of  Atchley  et  al.  (1990b).  Qualitatively,  the  numerical  results 
are  close  to  the  experimental  ones.  The  period,  1.95  ms,  is  identical  within  the  precision  with  which  it  can  be 
read  from  the  figure.  The  waveform  exhibits  a  strong  asymmetry,  with  the  negative  amplitude  much  smaller 
than  the  positive  one.  The  major  difference  between  calculations  and  experiment  is  the  amplitude,  which  is 
about  24.0  kPa  in  Fig.  2.10,  but  13.5  kPa  in  the  data.  In  order  to  match  to  measured  wave  amplitude  the 
temperature  difference  Th  —  Tc  should  be  decreased  by  about  25  K.  A  possible  explanation  for  this  difference 
may  be  the  following.  First  of  all,  the  temperature  values  reported  in  Atchley  et  al.  were  measured  at  the 
surface  of  the  tube,  rather  than  in  the  middle  of  the  heat  exchanger  plates.  For  the  conditions  of  the 
experiment  it  is  not  unreasonable  to  expect  a  temperature  difference  of  the  order  of  10-20  K  between  these 
two  points.  In  addition,  the  experimental  setup  most  likely  includes  several  loss  mechanisms  (e.g.  form  drag 
of  the  plates,  acoustic  streaming,  natural  convection)  not  included  in  our  model. 

2.6.2  Piston-Driven  Refrigerator 

Now  we  consider  a  simple  model  of  a  thermoacoustic  refrigerator  in  which  a  piston  at  the  left  end  of  the 
tube  sets  up  a  standing  wave;  the  right  end  is  modelled  as  rigid. 

The  driving  frequency  is  equal  to  the  natural  frequency  of  the  tube  open  at  one  end,  namely  cj  = 
27tVtRTV(4 £),  where  L  =  0.5  m  is  the  tube  length  and  I*  =  293  K  the  initial  uniform  temperature  of 
the  gas  and  solid  structures.  An  oscillating  velocity  u(t)  =  Ua  si  nut  is  prescribed  at  x  =  0  and  pl  =  p'  = 
V  =  0  there  as  discussed  at  the  end  of  section  2.2.  The  gas  is  helium  at  a  mean  pressure  P0  =307  kPa.  The 
stack  consists  of  0.28  mm-thick,  30  mm-long  plates  spaced  by  0.77  mm  with  a  thermal  conductivity  of  ks  = 
0.48  W/(K  m)  characteristic  of  fiberglass.  It  is  well  known  from  linear  theory  (see  e.g.  Swift  1988;  Karpov 
and  Prosperetti  1998)  that  the  thermoacoustic  heat  flux  is  strongest  when  the  stack  is  positioned  midway 
between  a  velocity  node  and  antinode.  Since  in  this  case  the  tube  length  is  one  quarter  of  the  acoustic 
wavelength,  we  position  the  stack  at  xs/L  =  0.5. 

The  hot  heat  exchanger  is  modelled  by  assuming  a  spatially  and  temporally  constant  temperature.  The 
cold  heat  exchanger  is  assumed  to  be  unloaded  and  to  only  exchange  heat  with  the  gas  and  the  stack.  It 
is  therefore  modelled  assuming  a  spatially  uniform,  time-dependent  temperature  the  mean  value  of  which 
will  decrease  with  time  under  the  action  of  the  sound  waves.  Since  the  thermal  conductivity  is  assumed 
to  be  very  large,  this  cold  heat  exchanger  is  simply  characterized  by  the  product  pscs.  The  duration  of 
the  transient  of  the  system  increases  with  the  value  of  this  quantity  and  therefore,  in  order  to  limit  the 
computational  time,  we  choose  pscs  =  480  kW/(I<  m3)  which  is  about  one  order  of  magnitude  smaller  than 
the  appropriate  value  for  realistic  materials.  In  spite  of  this  limitation,  the  results  that  follow  are  useful 
to  demonstrate  the  performance  of  the  model.  For  simplicity,  we  disregard  blockage  effects  and  take  both 
fv  and  fx  to  vanish  outside  the  stack  and  heat  exchanger  region  thus  neglecting  momentum  and  energy 
exchange  with  the  resonator  tube  walls. 

Figure  2.13  shows  the  temperature  Tw  as  a  function  of  time  at  the  cold  heat  exchanger  (lowest  line), 
at  the  stack  midpoint  (middle  line),  and  1  mm  away  from  the  end  of  the  hot  heat  exchanger.  Here  the 
imposed  velocity  amplitude  Ua  is  24.2  m/s  which,  when  converted  to  pressure  according  to  the  standard 
acoustic  relation  Ua  =  Pa/p^a ,  corresponds  to  a  drive  ratio  Pa/Po  —  0.04.  The  temperature  of  the  cold  heat 
exchanger  initially  decreases  with  time  and  finally  stabilizes  at  about  18.5  K  less  than  the  initial  temperature. 
At  steady  state,  on  average,  the  heat  extracted  from  the  heat  exchanger  by  the  gas  must  exactly  balance 
the  heat  gained  by  conduction  from  to  the  stack  plates.  Thus,  the  acoustic  power  supplied  to  the  system 
by  the  piston  is  spent  only  to  maintain  this  steady  state  by  removing  during  each  cycle  the  same  amount  of 
heat  from  the  cold  heat  exchanger  that  it  receives  from  the  stack  by  conduction.  Figure  2.14  shows  the  final 
mean  temperature  distribution  in  the  region  .of  the  stack  and  the  heat  exchangers.  The  mean  temperature 
of  the  solids  is  shown  by  the  solid  line,  and  that  of  the  gas  by  the  dotted  line.  There  are  near-jumps  in  the 
temperature  of  the  solid  structure  at  the  ends  of  the  stack  similar  to  those  found  in  Fig.  2.12.  Note  that 
a  significant  difference  between  the  gas  and  solid  structure  temperatures  occurs  only  in  the  heat  exchanger 
regions,  whereas  the  temperature  of  the  gas  and  the  plates  is  almost  the  same  over  most  of  the  stack.  The 
mean  temperature  of  the  gas  is  higher  than  that  of  the  solid  structure  over  one  part  of  the  heat  exchanger 
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Figure  2.13:  Mean  wall  temperature  as  a  function  of  time  at  the  cold  heat  exchanger  (dotted  line),  at  the 
stack  midpoint  (solid  line),  and  1  mm  away  from  the  hot  heat  exchanger  (dot-dash  line)  for  a  driven  tube 
with  a  drive  ratio  of  4%.  The  final  temperature  difference  between  the  heat  exchangers  is  18.5  K. 


Figure  2.14:  Mean  steady-state  temperature  distribution  in  the  region  of  the  stack  and  the  heat  exchangers 
for  the  solid  surfaces  (solid  line)  and  the  gas  (dotted  line)  for  the  driven  tube  of  the  previous  figure. 
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Figure  2.15:  Schematic  representation  of  a  thermoacoustic  prime  mover /thermoacoustic  refrigerator  combi¬ 
nation.  A  temperature  difference  T#  -  Tq  is  maintained  across  the  prime  mover  stack  by  its  heat  exchangers. 
The  standing  wave  generated  by  the  prime  mover  stack  causes  a  temperature  difference  across  the  refrigerator 
stack. 


and  it  is  lower  over  the  other  part.  The  same  result  was  found  in  the  two-dimensional  calculation  of  Worlikar 
(1997),  Worlikar  and  Knio  (1999),  and  Mozurkevich  (1998b). 

2.6.3  Thermoacoustic  Refrigerator  Coupled  with  Prime  Mover 

As  a  final  example  we  consider  the  combination  of  a  thermoacoustic  prime  mover  and  a  thermoacoustic 
refrigerator  sketched  in  Fig.  2.15  housed  in  a  1  m-long  rigidly  terminated  tube.  The  prime  mover  stack 
is  located  in  the  right  part  of  the  tube  at  xs/L  =  0.26  and  the  standing  wave  that  it  generate  induces  a 
temperature  difference  across  the  refrigerator  stack  positioned  in  the  left  part  of  the  tube  at  xs/L  ~  0.73. 

As  before  we  neglect  blockage  effects  and  gas-solid  momentum  and  energy  exchanges  away  from  the 
stack/heat  exchangers  region.  The  gas  is  helium  at  a  mean  pressure  of  307  kPa.  The  stacks  and  heat 
exchangers  consist  of  0.28  mm-thick  parallel  plates  spaced  by  0.77  mm  with  a  length  of  30  mm  and  1.5 
mm  respectively.  For  the  same  reasons  mentioned  above  we  set  the  product  pscs  to  480  kW/(K  m3).  The 
thermal  conductivity  of  the  stack  plates  is  ks  -  0.48  W/(K  m)  as  before,  whereas  the  plates  of  the  refrigerator 
heat  exchangers  are  assumed  to  have  infinite  thermal  conductivity.  The  temperature  of  the  prime-mover 
heat  exchangers  is  prescribed  to  be  Tc  =  293  K,  Tjj  —  493  K  and  held  fixed.  The  temperatures  of  the 
refrigerator  heat  exchangers  are  allowed  to  change  in  a  manner  similar  to  the  cold  heat  exchanger  of  the 
previous  example.  Initially  the  temperature  is  a  uniform  293  K  along  the  refrigerator  stack,  while  it  is  a 
linear  function  along  the  prime  mover  stack. 

Figure  2.16  shows  the  pressure  at  the  cold  end  of  the  tube  as  a  function  of  time.  There  is  an  initial  fairly 
rapid  rise  up  to  a  maximum  amplitude  reached  at  about  0.15  s,  followed  by  a  decline  to  the  steady  state 
regime  which  is  essentially  attained  at  0.4  s.  In  order  to  understand  this  non-monotonic  behavior  we  show 
in  Fig.  2.17  the  temperature  at  three  positions  in  the  prime  mover  stack,  the  midpoint  and  1  mm  away  from 
each  end.  We  see  here  that  the  temperatures  near  the  two  ends  of  the  stack  move  in  opposite  directions 
approaching  each  other  so  much  so  that  the  initial  temperature  gradient  of  6.67  K/mm  prevailing  along  most 
of  the  stack  is  reduced  at  steady  state  to  only  4  K/mm  (except  for  the  1  mm-long  segments  adjacent  the 
heat  exchangers).  The  amplitude  decrease  of  Fig.  2.16  clearly  occurs  as  a  consequence  of  this  trend.  While 
the  pressure  reaches  steady  state  at  about  0.4  s,  the  stack  temperatures  continue  to  slowly  adjust  after  this 
time  but  without  a  significant  change  in  the  temperature  gradient. 

Turning  now  to  the  refrigerator  section,  we  show  in  Fig.  2.18  the  temperatures  of  the  two  heat  exchangers 
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Figure  2.16:  Transient  behavior  of  the  pressure  at  the  cold  end  of  the  tube  for  the  thermoacoustic  prime 
mover/refrigerator  combination  sketched  in  the  previous  figure.  Note  the  non-monotonic  behavior  of  the 
pressure  amplitude. 


Figure  2.17:  Temperature  vs.  time  at  three  positions  in  the  prime  mover  stack  for  the  thermoacoustic  prime 
mover /refrigerator  combination  of  the  previous  two  figures.  The  dot-dash  and  dotted  lines  lines  are  for  points 
1  mm  away  from  the  hot  and  cold  heat  exchangers,  respectively;  the  solid  line  is  for  the  stack  midpoint.  The 
initial  temperature  gradient  of  6.67  K/mm  is  reduced  to  only  4  K/mm  at  steady  state. 
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Figure  2.18:  Temperature  vs.  time  of  the  cold  (dotted  line)  and  hot  (dot-dash  line)  heat  exchangers  and  the 
stack  midpoint  (solid  line)  in  the  refrigerator  unit  of  the  thermoacoustic  prime  mover/refrigerator  combina¬ 
tion  of  the  previous  three  figures. 


Figure  2.19:  Temperature  difference  between  the  two  heat  exchangers  vs.  time  in  the  refrigerator  unit  of  the 
thermoacoustic  prime  mover /refrigerator  combination  of  the  previous  figures. 
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Figure  2.20:  Mean  steady-state  temperature  distributions  in  the  the  refrigerator  stack  and  connected 
heat  exchangers  for  the  solid  surfaces  (solid  line)  and  the  gas  (dotted  line)  for  the  thermoacoustic  prime 
mover /refrigerator  combination  sketched  in  Fig.  2.15. 


Figure  2.21:  Mean  steady-state  temperature  distributions  in  the  the  prime  mover  stack  and  connected 
heat  exchangers  for  the  solid  surfaces  (solid  line)  and  the  gas  (dotted  line)  for  the  thermoacoustic  prime 
mover/refrigerator  combination  sketched  in  Fig.  2.15. 
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and  of  the  stack  midpoint,  all  as  functions  of  time.  The  temperatures  of  the  cold  and  hot  heat  exchangers 
initially  move  in  opposite  directions,  as  expected.  However,  after  about  1.0  s,  the  temperature  of  the  cold 
heat  exchanger  starts  increasing  due  to  viscous  and  thermal  heating  of  the  refrigerator  stack  and  its  heat 
exchangers.  For  the  refrigerator  stack  the  temperature  changes  much  more  slowly  than  for  the  prime  mover 
stack  reaching  its  steady  state  only  after  1.2  s.  Figure  2.19  shows  the  temperature  difference  between  the 
two  refrigerator  section  heat  exchangers  as  a  function  of  time. 

Figures  2.20  and  2.21  show  the  final  mean  temperature  distributions  for  the  refrigerator  and  prime  mover 
stacks  and  connected  heat  exchangers.  Significant  temperature  differences  between  the  plates  (solid  lines) 
and  the  gas  occur  only  near  the  heat  exchangers.  For  the  refrigerator  stack  the  gas  temperature  is  higher 
than  the  solid  temperature  over  one  part  of  the  heat  exchanger  and  lower  over  another  part  as  found  before. 

2.7  Conclusions 

We  have  described  a  simplified,  quasi-one-dimensional,  time-domain  mathematical  model  of  thermoacoustic 
prime  movers  and  heat  pumps.  The  model  is  quite  flexible  as  demonstrated  by  the  examples  that  we  have 
presented:  a  prime  mover,  a  combined  prime  mover /refrigerator  system,  and  a  piston-driven  refrigerator. 
We  have  been  able  to  follow  in  the  time  domain  the  evolution  of  these  systems  and,  for  the  first  two,  we 
have  presented  results  that  trace  the  evolution  of  the  initial  linear  instability  to  the  nonlinear  regime,  and 
finally  to  steady  state  where  the  instability  saturates  to  a  finite  amplitude.  This  is  a  specifically  nonlinear 
effect  that  cannot  possibly  be  captured  by  models  based  on  a  linear  approximation.  We  have  also  tested  the 
model  successfully  for  the  prime  mover/refrigerator  combination  assuming  periodicity  boundary  conditions; 
we  have  not  shown  these  results  as,  for  the  conditions  tested,  the  system  turned  out  to  be  stable. 

It  is  apparent  from  the  equations  presented  in  sections  2.2  to  2.4  that  the  model  includes  several  nonlinear 
mechanisms.  In  the  first  place,  it  accounts  for  mode-mode  coupling.  Even  though  in  actual  thermoacoustic 
engines  the  resonator  is  built  so  as  to  detune  the  harmonics,  at  sufficiently  large  amplitudes  a  coupling  still 
exists  and  represents  a  significant  mechanism  for  energy  loss:  in  a  thermoacoustic  engine,  typically  it  is  only 
the  fundamental  mode  that  is  unstable,  while  all  higher  modes  are  damped.  Hence,  energy  leaking  out  of  the 
fundamental  mode  into  higher  ones  is  dissipated;  a  weakly  nonlinear  analysis  of  this  phenomenon  is  presented 
in  Karpov  &  Prosperetti  (2000).  Similarly,  in  a  refrigerator,  the  stack  is  badly  positioned  with  respect  to 
the  spatial  distribution  of  pressure  and  velocity  for  modes  other  than  the  fundamental,  which  leads  to  a  loss 
of  efficiency.  Secondly,  the  model  does  not  contain  any  limitation  on  the  displacement  s  of  the  fluid  particles 
and,  therefore,  it  dispenses  with  the  standard  assumption  of  linear  theory  in  which  sdTw/dx  <C  Tw.  On  the 
other  hand,  it  must  be  recognized  that  other  nonlinearities  are  not  accounted  for.  For  example,  the  effects 
of  acoustic  streaming  or  those  of  vortex  shedding  from  geometric  discontinuities  are  not  included.  While  the 
latter  could  possibly  be  accounted  for  by  lumped  resistances  in  a  finite-difference  implementation,  inclusion 
of  the  former  appears  to  be  more  difficult.  It  may  also  be  noted  that  the  cross-stream  momentum  and  heat 
exchange  terms  are  modelled  linearly,  which  is  however  justified  by  the  fact  that,  typically,  the  stack  is  many 
hydraulic  diameters  long,  so  that  a  fully  developed,  parallel-flow  approximation  is  warranted. 

While  these  and  other  mechanisms  and  other  details  of  a  real  thermoacoustic  system  may  not  be  easily 
incorporated  in  the  model,  a  great  variety  of  design  parameters  can  be  accounted  for  such  as  geometry 
(and  in  particular  a  non-constant  cross  section),  physical  properties,  temperature  conditions,  standing  and 
travelling  waves,  and  others.  Thus  we  believe  that  -  at  the  very  least  -  the  present  model  can  be  used  in  a 
relative  sense  to  compare  the  effect  of  proposed  design  modifications  and  to  gain  insight  into  the  performance 
of  new  systems. 
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Chapter  3 

Linear  and  Weakly  Nonlinear  Theory 
of  Thermoacoustic  Devices 


The  nonlinear  model  described  in  the  preceding  chapter  and  its  numerical  treatment  are,  we  believe,  the 
most  significant  piece  of  work  in  Thermoacoustics  completed  in  the  course  of  this  project.  In  parallel,  we 
have  also  conducted  analytical  work  aimed  at  finding  a  time-domain  description  of  the  phenomena,  both 
in  the  linear  and  in  the  nonlinear  regime.  This  work  has  been  based  on  the  same  approximate  formulation 
described  in  the  previous  chapter  but,  since  the  different  normal  modes  are  here  treated  individually,  the 
terms  responsible  for  the  cross-axis  exchange  of  momentum  and  energy  can  be  treated  more  accurately. 

Here  we  only  present  a  brief  description  of  this  work;  full  details  are  given  in  the  two  publications  resulting 
from  these  studies:  Karpov  &  Prosperetti  (1998)  and  Karpov  &  Prosperetti  2000. 


3.1  Introduction 


This  work  is  based  on  a  perturbative  treatment  of  the  problem,  which  Is  motivated  by  the  observation  that 
both  experiment  (see  e.g.  Wheatley  1986;  Atchley  1994)  and  the  calculations  described  in  the  previous 
chapter  show  that  the  initial  build-up  of  the  thermoacoustic  instability  has  the  character  of  a  modulated 
standing  wave  the  frequency  of  which  is  essentially  dictated  by  the  resonator,  while  the  amplitude  is  slowly 
varying  in  time.  Similarly,  if  the  thermoacoustic  device  is  driven  externally  by  a  loudspeaker  or  a  piston, 
the  steady  state  temperature  distribution  evolves  slowly  over  the  time  scale  of  the  acoustic  period.  These 
observations  suggests  the  possibility  of  setting  up  a  perturbation  scheme  based  on  the  smallness  of  the  ratio 
of  the  characteristic  period  of  oscillation  to  the  characteristic  time  for  the  evolution  of  the  thermoacoustic 
effect. 

We  start  from  the  same  mathematical  model  of  the  previous  chapter,  which  is  here  repeated  for  the 
reader’s  convenience: 

-o. 


di  S  dx 
dpu  1  d  /ri 

ir + 


dp 

dx 


=  - T>(pu ), 


(3.1) 


(3.2) 


pcpn  Tw-T)-cp 


(3.3) 


All  the  terms  are  described  in  detail  in  section  2  of  the  previous  chapter;  time  is  denoted  by  t  for  reasons  that 
will  become  soon  apparent.  For  simplicity,  and  in  view  of  the  large  thermal  capacity  of  typical  stack  plates, 
we  assume  in  this  chapter  that  Tw  =  Tw(x)  is  a  prescribed  function  of  x.  With  relatively  straightforward 


1-37 


adjustments,  however,  the  present  method  can  be  extended  to  deal  with  a  Tw  dependent  on  time  over  the 
fast  time  scale  of  the  acoustic  oscillations  as  well  as  over  the  slower  time  scale  of  the  developing  instability. 

For  the  exchange  terms  V,  and  Q  we  use  the  exact  form  of  the  linear  theory  as  described  in  the 
previous  chapter,  namely 

V{pu)  =  itu -  — —  pou ,  (3.4) 

1  -  Iv 

H(TW  -T)  =  (Tw  ~  T) ,  (3.5) 

1  ~  JK 


Q(pu)  =  ^ 


fv 


1  -fv 


<rfi<  \ 
1  -fid 


PoU, 


(3.6) 


where  the  functions  fvj<  depend  on  the  viscous  and  thermal  penetration  length  at  the  frequency  cj  and  are 
given  in  (2.22)  for  the  plane  geometry  and  in  (2.30)  for  the  cylindrical  geometry.  It  will  be  noted  that  we 
have  replaced  the  density  p  iun  the  right-hand  side  by  the  unperturbed  equilibrium  density  po*  This  we  do 
for  simplicity  on  the  basis  of  numerical  results  obtained  by  the  method  of  the  previous  chapter  which  show 
that  this  simplification  has  a  negligible  effect. 

In  the  model  (3.1)  to  (3.3),  the  terms  in  the  left-hand  sides  describe  non-linear  oscillations  in  a  gas 
column  with  variable  cross  sectional  area  and  temperature  stratification  and  are  therefore  responsible  for 
the  “carrier”  frequency  of  the  wave.  The  heart  of  the  thermoacoustic  effect  is  in  the  terms  in  the  right-hand 
sides.  The  observed  slowness  of  the  modulation  implies  that  the  effect  of  these  terms  over  a  single  period  of 
oscillation  is  small.  As  will  be  seen  in  the  following,  this  effect  arises  through  integrals  over  the  tube  length. 
Therefore,  the  effect  will  be  small  not  only  when  the  terms  themselves  are  small  but  -  as  usually  happens  in 
practice  -  when  they  are  large  only  over  a  small  fraction  of  the  tube  length.  This  statement  can  be  formally 
justified  as  shown  in  Appendix  A  of  Karpov  &  Prosperetti  (2000). 

In  order  to  set  up  a  perturbation  scheme,  we  introduce  a  small  parameter  e  and  set 


Fk  = 


fv  ~  f_K 

(l-/v)(l-£r)  ’ 


fv 

1  -fv’ 


(3.7) 


with  Fk,q  formally  treated  as  0(1)  quantities.  As  shown  in  Karpov  &  Prosperetti  (2000),  this  does  not 
necessarily  imply  that  the  /’s  are  of  order  e,  but  only  that 

-j-  fv,i<  ~  €,  (3.8) 

where  Ls  is  the  length  of  the  stack  region,  i.e.,  the  region  where  the  /’s  are  not  small.  An  explicit  definition 
of  6  is  not  necessary  as  the  final  results  do  not  explicitly  depend  on  this  parameter,  but  one  may  think  of  it 
as  the  ratio  of  the  standing  wave  period  to  the  modulation  time  scale.  We  use  the  method  of  multiple  time 
scales  (see  e.g.  Kevorkian  and  Cole  1996;  Murdock  1991;  Hinch  1991)  and  introduce  the  new  time  variables 


t  —  t,  r  —  ei, 


0  =  e2£,  rj  —  e3t . 


(3.9) 


As  a  consequence  of  these  definition  we  have 


j9  _  d_ 

di  dt+€  dr 


4-  e 


89 


3  9 

+  e3— +  ••• 


drj 


(3.10) 


The  field  variables  are  also  expanded  in  a  power  series  in  e;  for  example 


p{x.t)  -  po  +  epi(x,  t)  +  e2p2{x,  i)  4-  e3p3(x,  t)  4-  e4p4(x,  £)  +  ...  ,  (3-11) 


with  analogous  expressions  for  v!  etc.  These  expansions  imply  that  the  nonlinearity  is  taken  to  be  of  the 
same  order  as  the  amplitude  modulation,  which  is  the  interesting  case.  Indeed,  if  the  modulation  were 
much  stronger  than  the  nonlinearity,  we  would  essentially  be  dealing  with  the  linear  problem  (see  Karpov  & 
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Prosperetti  1998  for  an  analysis  of  it) .  On  the  other  hand,  since  the  wave  exists  only  due  to  the  thermoacoustic 
instability,  one  cannot  have  a  strong  nonlinearity  coupled  with  weak  thermoacoustic  effects.  As  a  consequence 
of  an  expansion  analogous  to  (3.11)  for  u,  we  may  write 

V{pu)  -  e2T>i  4-  esV2  +  e4X>3  +  -  - .  ;  (3.12) 


the  expansion  begins  with  a  term  of  order  e2  since  V  is  itself  of  order  e  as  implied  by  (3.4)  and  (3.7). 

For  clarity,  it  is  important  to  stress  two  important  aspects  in  which  the  work  described  in  this  paper 
differs  from  most  other  non-linear  stability  studies.  In  the  first  place,  while  the  deviation  from  marginal 
stability  conditions  is  usually  measured  by  a  single  control  variable,  here  it  is  a  whole  function  -  the  wall 
temperature  distribution  Tw  -  that  determines  the  stability  properties  of  the  linear  system.  The  situation 
might  be  reduced  to  the  more  usual  one  by  assuming  that  the  temperature  distribution  in  the  stack  has  a 
certain  functional  form  dependent  upon  one  parameter  which  would  then  play  the  role  of  control  variable. 
For  example,  a  linear  temperature  distribution  would  be  characterized  by  the  temperature  gradient.  Such 
an  assumption  is  however  unnecessarily  restrictive  and  it  is  preferable  to  keep  the  framework  general. 

Secondly,  perturbation  expansions  are  usually  carried  out  in  the  neighborhood  of  linear  marginal  stability 
conditions  that  are  known  exactly.  In  order  to  proceed  in  this  way,  we  would  have  to  solve  the  linear  problem 
including  the  exchange  terms  fv,K,  find  the  marginal  stability  conditions,  and  then  allow  for  a  perturbation. 
Here  such  a  procedure  would  not  lead  to  very  transparent  results  given  the  complexity  of  the  linear  problem, 
compounded  by  the  fact  that,  with  non-zero  /’ s,  the  linear  operator  is  not  self-adjoint.  Hence,  we  prefer  to 
take  a  different  approach  choosing  as  our  base  state  the  situation  in  which  there  is  no  coupling  between  the 
gas  and  the  solid  structure  even  though,  in  general,  this  situation  is  not  marginally  stable.  The  marginal 
stability  condition  will  be  determined  as  part  of  the  perturbation  procedure  itself  much  as  in  our  earlier  paper 
(Karpov  and  Prosperetti  1998)  of  which  the  first  few  steps  of  the  present  perturbation  procedure  constitute 
a  refinement. 

Here,  therefore,  as  the  parameter  e  is  increased  from  zero,  two  effects  appear:  the  coupling  with  the  solid 
structure  dampens  or  amplifies  the  oscillations,  and  nonlinear  effects  influence  their  development.  It  will 
be  seen  that  the  first  effect  sets  in  at  a  lower  order  in  e  than  the  second  one.  The  first  few  steps  of  the 
perturbation  procedure,  therefore,  furnish  successively  closer  approximations  to  the  linear  stability  threshold, 
from  which  the  nonlinear  effects  eventually  take  off. 

In  order  to  deal  with  this  situation  we  must  clarify  the  role  played  by  the  wall  temperature  distribution 
Tw(x).  In  the  usual  perturbation  procedures,  the  small  parameter  is  the  difference  between  the  given  value 
of  the  control  variable  and  the  critical  value.  Here  we  deal  with  the  function  TW}  and  the  critical  conditions 
are  not  known  exactly.  Thus  we  represent  the  given  function  Tw  as 

Tw(x)  =  Twq(x)  -f-  eTwi(x)  +  €~TW2(x)  4- . . .  .  (3.13) 

The  terms  Tw 0,  Tw\  will  be  determined  so  as  to  satisfy  the  marginal  stability  conditions,  after  which  the 
difference  Tw(x)  —  [Two(x)  4-  eTw\(x)]  will  be  the  driving  force  for  the  development  of  the  instability  to  the 
order  considered.1  In  principle,  the  same  procedure  can  be  applied  to  higher  orders  in  e.  Since  the  appearance 
of  the  instability  depends  on  integral  conditions  on  the  Twj’ s,  in  principle,  there  is  a  degree  of  non-uniqueness 
here  since  the  Twj  ’s  can  be  taken  arbitrarily,  provided  only  the  integral  conditions  are  satisfied.  In  spite  of 
this  seeming  level  of  arbitrariness,  one  would  expect  the  results  to  be  insensitive  to  any  specific  choice,  and 
indeed  we  have  verified  this  numerically  as  will  be  described  below. 

From  (3.13),  we  also  have 

~ -  =~Go  +  cGi  +  e2G2  +  . . .  .  (3.14) 

ax 

where 

=  J  =  0,1,2,....  (3.15) 

xMore  precisely,  one  could  determine  a  value  of  T*o2  corresponding  to  marginal  stability,  after  which  the  difference  Tw(x)  — 
+  eTwi(x)  +  e2T£2]  would  be  the  driving  force  for  the  instability. 
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Furthermore,  since  the  equilibrium  density  is  connected  to  p0  and  Tw  by  the  equation  of  state,  when  we 
write 

p(x,t)  =  po  +  epi(x,£)  +  e2p2{x,t)  +  esp3{x,i)  +e*p4(x,i)  +  ...  ,  (3.16) 

the  corrections  pi  etc.  will  contain  time-independent  components  related  to  Tw i  etc.  It  is  also  convenient  to 
introduce  a  term 

T(p,  u)  =  pH  (Tw  -T)-^  Q{PU)  =  -iufK(p  -  Po)  -  pcp  (1  ’  (3-17) 

where  the  last  step  is  justified  by  relations  valid  in  the  linear  approximation.2  Upon  substituting  the  previous 
e-expansions,  we  find 

T{p ,  u)  =  e2Ti  +  e3T2  +  e4T3  +  . . .  ,  (3.18) 

where 

T\  =  -iojFkPi  —poCpFgGoUi ,  V i  —  ip0ojFyUi ,  (3.19) 


etc. 

All  the  expansions,  together  with  (3.10),  are  now  substituted  into  the  equations  (2.1)  to  (2.11)  and  the 
various  orders  in  e  are  separated  giving  rise  to  a  sequence  of  problems.  At  order  zero  the  solution  is  simply 

«o  =  0,  (3.20) 

with  po  and  T0  =  Two{x)  arbitrary  but  regarded  as  given. 


3.2  The  First-Order 

At  order  e  the  equations  are 


Problem 


dpi  1  d 

-m+saiiSpoUl)  =  0’ 


dui  dpi 

Po~dT  +  ~d^  ~  °’ 


dpi  IPo  d(Sui)  _ 

dt  s  dx 


from  which  the  pressure  field  is  found  to  be 

Pi  =  A(r,  0,  p)  Pi  ( x )  exp  ( iut )  +  c.c. , 


(3.21) 

(3.22) 

(3.23) 

(3.24) 


where  A  is  the  slowly-varying  amplitude,  c.c.  denotes  the  complex  conjugate,  and  Pi  is  the  solution  of 


(3.25) 


where 

co(x)  =  7  FTw0(x),  (3.26) 

is  the  local  adiabatic  sound  speed.  We  take  the  tube  to  be  rigidly  terminated  at  the  ends  and  therefore 
impose  that 

=  0  at  x  =  0,  x  =  L ,  (3.27) 

dx 

2We  have  verified  numerically  that  the  use  of  this  relation  in  the  nonlinear  regime  leads  to  a  better  agreement  with  the 
complete  nonlinear  theory. 
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which  ensure  that  the  eigenvalue  co2  is  real  and  positive  (see  e.g.  Morse  and  Feshbach  1953,  p.  728;  Naylor 
and  Sell  1982,  p.  502).  The  eigenfunction  Pi  can  therefore  also  be  taken  real  and,  for  later  convenience,  we 
normalize  it  so  that 

[  S(x)P2(x)  dx=  Vp2,  (3.28) 

Jo 

where  V  is  the  volume  of  the  device. 

The  way  in  which  Tw 0  is  to  be  understood  needs  some  clarifications,  that  amplify  the  comments  made 
before  Eq.  (3.13).  Let  us  consider  first  the  case  in  which  a  certain  temperature  distribution  Tw(x)  is 
prescribed.  If  one  plans  to  study  only  the  linear  problem,  Tw0  can  be  taken  as  the  given  Tw.  On  the  other 
hand,  if  the  plan  is  to  carry  the  expansion  to  include  nonlinear  effects,  as  will  be  discussed  below,  Tw0  must 
be  determined  (or,  at  least,  constrained)  using  the  results  of  the  next  step  in  the  perturbation  procedure. 
Finally,  if  the  objective  is  the  determination  of  the  onset  temperature  distribution,  Tw(x)  is  unknown  at  the 
outset,  and  the  Tw 0  appearing  in  (3.25)  is  its  first  approximation  that  will  be  determined  at  the  next  step. 
Actually,  the  last  two  possibilities  are  one  and  the  same  as,  if  one  plans  to  carry  the  expansion  to  include 
nonlinear  effects,  one  needs  to  determine  the  onset  temperature  distribution  first.  A  consequence  of  the 
indeterminacy  of  Tw0  at  this  stage  is  an  indeterminacy  of  p0,  which  will  be  determined  from  the  equation  of 
state  7ZpoTwo  =  Po  once  Tw0  is  found. 

In  deriving  an  expression  for  the  gas  density  to  this  order,  we  anticipate  the  fact  that,  at  the  next  order, 
we  will  encounter  a  contribution  Twl  that  will  affect  the  undisturbed  value  of  p.  Hence  we  write 

Pi  =  pio  +  Pn  ,  (3.29) 


where  pio  will  be  determined  at  the  third  step  from  the  equation  of  state 

1Z  (po  +  epio)  (Tw o  +  tTwi)  =  po  . 


(3.30) 


The  term  pu  is  instead  associated  to  pi  and  can  be  found  from  (3.21)  to  (3.23),  together  with  m.  The 
dependence  of  both  quantities  on  the  time  variables  and  the  perturbation  amplitude  A  is  the  same  as  (3.24), 
and  the  spatial  dependence  is  given  by 


i?i(x)  =  - 


1  d 
l d'2S  dx 


(3.31) 


Ui(x)  - 


i  dP\ 
to  po  dx 


(3.32) 


3.3  Second-Order  Problem 

At  order  e2,  one  finds  the  following  equation  for  po: 


d2P2  _  ld_  (  .2  qdP2 
dt 2  S  dx  I  0  dx 


where  Cq  is  defined  in  (3.26)  and 


frfrro x  _  /  .x  g7l  9 

(rhs)2  (/  i)  dt  2  9tdr 


d_ 

dt 


=  (. rhs)2 , 


dpi  ,  7Pi  9 


^  7Po  9  - 

+  S  dxb 
fn.  9A  A 

=  ~ 


1  -r,  ,  dui  pi  dui 

—  T>i+  «i-q — I - jr—  , 

Po  ox  po  dt  J 

1  A 

S  dx 
fv  -  fl< 


^L-ujpoSc?ouPj  +  (7  -  1  )u2fKPx 


(opoUi  dTw o 


(1  —  <r)(  1  —  fv)  1  Tw o  dx  J 


exp  iwt  +  c.c. . 


(3.33) 


(3.34) 
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The  forcing  at  frequencies  ±u  in  the  right-hand  side  will  generate  resonant  terms  proportional  to 
t  exp  (±iojt)  in  the  solution  for  p2  which  would  lead  to  a  breakdown  of  the  approximation  over  times  of 
order  (ecu)-1.  To  avoid  this  resonance,  as  in  the  standard  procedure  (Kevorkian  and  Cole  1996;  Murdock 
1991;  Hinch  1991),  we  impose  the  solvability  condition  that  the  right-hand  side  of  the  equation  be  orthogonal 
to  the  solutions  of  the  (adjoint)  homogeneous  equation,  namely  exp  (±iuit)  Pi : 


P'Ztt/uj  rl 

/  d*  / 

Jo  Jo 


S  [exp(=pzo;t)Pi]  (RHS)2  dx  —  0. 


(3.35) 


Multiplication  by  S(x)  before  integration  is  necessary  so  that  the  x-operator  in  the  left-hand  side  of  (3.33) 
be  self-adjoint  with  the  boundary  conditions  (3.27).  The  two  conditions  (3.35)  give 


8A 

dr 


i£l\  .4  =  0, 

OT 

plus  its  complex  conjugate,  where 

+  (7  -  1  )u2FkPx  +  (7  -  1)cvFqG0P\ 


dPi 

dx 


dx . 


(3.36) 


(3.37) 


The  solution  of  (3.36)  is 

A(r,ff,rj)  =  5(0,77)  exp  (iftir)  ,  (3-38) 

which  increases  exponentially  if  ImHi  <  0.  It  will  be  noted  that  Imfti  exists  only  due  to  the  P-terms  that 
account  for  momentum  and  energy  exchange  with  the  solid  structure  and  is  therefore  not  affected  by  the 
indeterminacy  of  pio  at  this  stage  of  the  calculation.  It  will  be  seen  below  that,  in  order  to  continue  the 
perturbation  scheme  into  the  nonlinear  domain,  for  consistency  it  is  necessary  to  have  Im  Cli  =  0,  which  is  a 
constraint  on  G0  and,  therefore,  Tw 0.  This  constraint  can  of  course  be  met  in  an  infinity  of  ways  although  one 
would  not  expect  any  particular  choice  to  have  consequences  for  the  remainder  of  the  calculation  provided 
the  choice  is  such  that  Tw  -  Tw0  =  0(e)  in  keeping  with  (3.13).  Indeed,  as  will  be  shown  below,  this  is  in 
agreement  with  our  numerical  evidence. 

Although  essentially  a  linear  result,  (3.37)  is  already  rather  interesting.  We  have  shown  (Karpov  & 
Prosperetti  1998)  how,  assuming  a  linear  temperature  gradient,  one  can  deduce  from  it  the  onset  value  Go- 
By  means  of  the  short-stack  approximation,  it  was  also  possible  to  derive  simple  formulae  that  generalize 
the  concept  of  critical  gradient  of  elementary  theory  to  the  presence  of  viscous  effects  and  narrow  gaps. 

When  (3.35)  is  satisfied,  one  can  look  for  the  solution  of  Eq.  (3.33)  for  p2  in  the  form  of  the  superposition 
of  terms  at  frequencies  0,  ±u),  and  ±2lo  respectively: 

P2  =  AP21(x)  eiu)t  +  A2P22(x)  e2iut  +  A*AP20(x)  +  c.c. .  (3.39) 

Here,  the  asterisk  denotes  the  complex  conjugate.  The  equations  for  P2 i(x),  P22(x ),  and  P2o(x)  are  obtained 
from  Eq.  (3.33)  upon  separating  the  different  frequency  components  in  the  right-hand  side.  However  P2o  is 
more  easily  determined  from  the  momentum  equation  (2.53).  The  explicit  form  of  these  equations  is  given 
in  Eqs.  (B7)  to  (B9)  in  Appendix  B  of  Karpov  &  Prosperetti  (2000). 


3.4  The  Landau  Equation 


The  procedure  demonstrated  in  the  previous  section  can  be  continued  to  order  e4  to  find  an  evolution 
equation  for  the  amplitude  A  of  the  well-known  Landau  form.  We  write  this  equation  in  terms  of  the 
unsealed  amplitude  A  =  e.4  exp(iut): 


dA 

dt 


j^ca  +  "I"  +  ^4-2  A  ^3)  \A\" 


A, 


(3.40) 
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5ji  =0.34 


Figure  3.1:  Negative  of  the  instability  growth  rate,  i.e.  imaginary  part  of  the  exact  linear  eigenfrequency  (solid  line), 
compared  with  the  successive  approximations  Imf^i  (dotted  line),  Im(fii  -i-fils)  (short  dashes),  and  Im(Hi  +  ^3) 
(long  dashes)  as  functions  of  the  normalized  position  xs/L  of  the  stack  for  the  geometry  described  in  the  text. 


where  is  obtained  by  multiplying  (3.37)  by  e, 


Hi 


1 

2V PqW 


fv  ,  Pio  \ 
1  -  fv  Po  ) 


+  (7  -  1)  w2  fl<  P\  + 


(7  -  1  )(fv  -  M 
(1  -  a)(  1  -  fv) 


cp  G  0  Pi 


dPi 

dx 


(3.41) 


with  similar  expressions  for  Ojt  =  e*fl/ t;  the  A’s  are  also  integrals  over  the  volume  of  the  thermoacoustic 
device  and  are  given  explicitly  in  Karpov  &  Prosperetti  (2000). 

Equation  (3.40)  contains  both  linear  and  nonlinear  contributions  that  it  is  better  to  discuss  separately. 
In  all  the  numerical  examples  that  follow  we  use  as  reference  case  the  geometry  of  the  experiment  of  Atchley 
et  al.  (1990),  the  implementation  of  which  in  the  context  of  the  present  mathematical  model  is  described  in 
detail  in  an  earlier  paper  (Yuan  et  al.  1997).  Briefly,  the  system  consists  of  a  38.2  mm-diameter  tube  with 
a  length  of  99.89  cm,  a  35  mm-long  stack  of  stainless  steel  plates  located  87.95  cm  from  the  cold  end,  and 
two  heat  exchangers.  The  operating  pressure  was  po  —  307  lcPa. 


3.4.1  Linear  regime 

In  the  linear  approximation,  Eq.  (3.40)  is 

~  =  i  (w  +  fh  +  Cl2  +  63)  A  ■  (3-42) 

In  addition  to  giving  the  growth  rate  of  the  perturbation  for  a  given  temperature  distribution  Tw,  Eq.  (3.42) 
can  be  used  to  determine  onset  conditions  by  adjusting  the  temperature  distribution  in  such  a  way  that  Im 
(&i  t  n*2  ~h  1^3)  =  0. 
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We  have  shown  (Karpov  &  Prosperetti  1998)  that,  retaining  just  the  first  term  the  critical  temperature 
gradient  found  in  this  way  is  already  in  excellent  agreement  with  Rott/s  exact  linear  theory  prediction.  A 
typical  example  is  shown  in  Fig.  3.1,  where  the  imaginary  part  of  the  exact  linear  eigenfrequency  -  i.e. 
minus  the  growth  rate  of  the  instability  -  (solid  line)  is  compared  with  the  successive  approximations  Im 
Hi  (dotted),  Im  (Oi  +  fi2)  (dashed),  and  Im  (fii  +  fi2  +  ^3)  ( long  dashes).  Here  the  horizontal  axis  is  the 
normalized  position  xs/L  of  the  cold  end  of  the  stack.  Further  results  of  this  type  will  be  found  in  Karpov 
&  Prosperetti  (1998). 

Results  such  as  these  confirm  the  correctness  of  the  perturbation  technique  and  the  accuracy  of  the 
analysis. 

We  now  turn  to  the  more  interesting  nonlinear  case. 

3.4.2  Nonlinear  regime 

Equation  (3.40),  and  its  complex  conjugate,  predict  the  evolution  of  the  oscillation  amplitude  in  time.  These 
equations  can  be  combined  to  find 


d|J-  =  -2  [lm03  -f  ImA3|A|2J  \A\2  ,  (3.43) 

which  is  readily  solved  with  the  result 

|I|2  =  - r . 1% - r - r  ,  (3-44) 

1  +  exp  (2 n3ii)  (\a\U\a\2o  - 1) 


where  Aq  denotes  the  initial  value  of  the  amplitude  and 


\A\lat 


im  n3 
Im  A3 


i 

A  3f 


(3.45) 


Since  f \zi  <  0  for  an  unstable  system,  this  equation  shows  that,  for  t  -¥  00,  \A\ 2  is  asymptotic  to  the 
saturation  value  |A|^,  as  also  follows  directlyfrom  the  right-hand  side  of  (3.40)  since  only  fi3  and  A3  have 
non-zero  imaginary  parts.  In  the  stable  case,  ft3;  is  positive  and  the  amplitude  decays  to  zero. 

Figure  3.2  shows  the  time  dependence  of  the  cold-end  pressure  for  the  middle- temperature  case  AT  = 
346  K  of  Atchley  et  al.  (1990)  and  exhibits  the  typical  saturation  of  the  instability  due  to  non-linear  effects. 

The  result  (3.45)  for  the  asymptotic  saturation  amplitude  is  of  particular  interest  as  it  contains  the  effect 
of  a  large  number  of  variables  such  as  the  deviation  of  the  temperature  distribution  from  onset  conditions, 
the  shape  of  the  resonator,  the  effect  of  stack  geometry  (through  the  values  of  the  exchange  parameters  fvj<) 
and  stack  length,  Prandtl  number,  and  others. 

We  assume  a  linear  temperature  variation  between  a  cold  temperature  Tc  and  a  hot  temperature  Th  in 
the  stack,  and  constant  temperatures  equal  to  Tc  and  Th  at  the  left  and  right  of  the  stack,  respectively. 

The  solid  line  in  Fig.  2  is  the  (positive)  peak  amplitude  as  a  function  of  the  temperature  difference  along 
the  stack  for  the  experimental  situation  of  Atchley  et  al.  (1990).  The  dotted  line  shows  |Asat|*3  For  this  case 
our  calculations  indicate  an  onset  temperature  of  573.5  K  which,  with  Tc  —  293  K,  corresponds  to  an  onset 
temperature  gradient  of  8.01  K/mrn.  One  recognizes  here  the  typical  structure  of  a  supercritical  bifurcation. 
In  this  figure  the  symbols  represent  the  three  data  points  reported  by  Atchley  et  al.  (1990).  In  order  to 
reconcile  these  data  with  our  theoretical  prediction,  one  would  have  to  assume  hot-end  temperatures  36,  22, 
and  20  K  lower  than  the  reported  values.  We  found  a  similar  effect  in  our  earlier  numerical  study  (Yuan 
et  al.  1997).  Two  possible  concurrent  explanations  come  to  mind.  The  first  one  is  the  temperature  jump 
that  is  established  between  the  ends  of  the  stack  and  the  heat  exchangers,  which  has  the  effect  of  reducing 
the  temperature  gradient  in  the  stack  (Brewster  et  al.  1997).  Secondly,  the  temperature  values  reported  in 

3Since  the  pressure  perturbation  is  the  sum  of  two  complex  conjugate  quantities,  the  dotted  line  actually  shows  2  [Re,4sa#|- 
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(a)  AT  -  346  K 


time  (sec) 


(b)  AT  =  346  K 


time  (sec) 


Figure  3.2:  Time  dependence  of  the  cold-end  pressure  for  the  middle-temperature  case  AT  =  346  K  of 
Atchley  et  al.  (1990b)  as  predicted  by  Eq.  (3.11)  with  each  term  computed  using  the  A  given  by  (3.40). 
The  lower  figure  is  an  enlargement  of  the  initial  phase  of  the  upper  one;  the  dashed  line  is  the  envelope  from 


Figure  3.3:  The  solid  line  shows  the  peak  amplitude  of  the  pressure  perturbation  at  the  cold  end  of  the 
resonant  tube  as  a  function  of  the  temperature  difference  AT  along  the  stack;  the  dashed  line  is  the  saturation 
amplitude  of  the  fundamental  mode  \A\sat  from  Eq.  (3.45);  the  symbols  are  the  data  points  reported  by 
Atchley  et  al.  (1990b). 


Atchley  et  al.  were  measured  at  the  surface  of  the  tube,  rather  than  in  the  stack.  As  mentioned  in  Yuan  et 
al.  (1997),  for  the  conditions  of  the  experiment  one  may  reasonably  expect  a  temperature  difference  of  this 
order  between  these  two  positions.  In  addition,  the  experimental  setup  will  most  likely  include  several  losses 
(e.g.  due  to  the  lack  of  alignment  between  the  heat  exchanger  and  stack  plates,  which  causes  an  increase 
of  the  effective  blockage  of  the  tube,  heat  conduction  in  the  gas,  and  others)  that  are  not  included  in  our 
model.  Hence  the  difference  with  the  data  is  in  the  expected  direction  and  of  a  reasonable  magnitude.  The 
discrepancy  at  the  lowest  temperature  difference  is  somewhat  larger  but,  on  the  other  hand,  it  is  here  that 
the  effect  of  losses  would  be  greatest. 

The  expressions  derived  before  contain  a  dependence  on  many  variables  the  effect  of  which  can  also  be 
studied.  Several  examples  are  given  in  Karpov  &  Prosperetti  (2000). 


3.5  Conclusions 

We  have  presented  a  time-dependent  weakly  nonlinear  theory  of  the  build-up  of  unstable  oscillations  in  a 
thermoacoustic  prime  mover.  The  expressions  we  have  derived  simulate  the  initial  growth  of  the  instability 
as  well  as  its  subsequent  saturation.  The  results  are  in  reasonable  agreement  with  the  limited  number  of 
experimental  observations  reported  in  the  literature  for  the  configuration  that  we  can  model. 

For  all  the  examples  that  we  have  discussed,  a  solution  of  Rott’s  equation  shows  that  the  imaginary 
part  of  all  the  modes  higher  than  the  fundamental  one  is  positive,  so  that  they  are  all  stable.  Therefore, 
the  second  and  higher  harmonic  amplitudes  are  non-zero  only  due  to  non-linear  energy  transfers  from  the 
fundamental.  In  the  present  perturbation  scheme,  energy  transfer  to  the  second  mode  starts  occurring  at 
the  second  order.  However,  the  energy  loss  of  this  mode  due  to  the  coupling  with  the  solid  structure  only 
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arises  at  the  fourth  order.  Since  the  amplitude  of  the  fundamental  saturates  at  a  level  such  that  the  energy- 
transferred  to  the  second  mode  can  be  dissipated  by  the  mechanisms  affecting  this  mode,  in  order  to  find 
the  saturation  amplitude  of  the  fundamental,  we  had  to  carry  the  expansion  to  this  order.  The  fact  that  the 
perturbation  method  only  includes  a  small  number  of  modes  prevents  energy  transfer  from  the  fundamental 
to  the  higher  modes.  As  a  consequence,  our  results  tend  to  overestimate  somewhat  the  saturation  amplitude. 
Nevertheless,  the  comparison  of  section  5  shows  that  there  is  a  parameter  region  where  the  error  is  acceptable. 

Even  if  our  results  may  not  be  sufficiently  precise  in  absolute  terms  when  the  temperature  gradient  is 
too  far  from  the  onset  value,  they  are  still  useful  for  comparative  purposes  as  they  indicate  trends  of  the 
system  performance  when  design  parameters  or  operating  conditions  are  varied.  For  example,  we  find  a  very 
strong  effect  of  cross-sectional  area  variations.  Since  the  expressions  we  have  derived  contain  a  great  amount 
of  detail  on  the  geometry  and  other  characteristics  of  the  system,  similar  sensitivity  studies  with  respect  to 
other  parameter  variations  are  possible. 

In  spite  of  the  algebraic  complexity  of  the  final  form  of  the  results,  the  present  solution  is  certainly  much 
easier  to  evaluate  than  carrying  out  a  full-fledged  multi-dimensional  ab  initio  numerical  calculation.  In  our 
earlier  paper  (Karpov  &  Prosperetti  1998)  it  was  possible  to  simplify  the  form  of  the  solution  by  using 
the  short  stack  approximation.  Here  this  approximation  is  not  very  useful  because  the  result  depends  on 
auxiliary  fields  such  as  P21 ,  P22 ,  etc.  the  determination  of  which  requires  the  solution  of  equations  which  can 
only  be  solved  in  special  cases.  Nevertheless  it  might  be  possible  to  develop  some  comparable  approximation. 
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Chapter  4 

Bubble  Screens 


During  the  last  few  months  of  the  ONR  grant,  we  have  also  worked  on  a  subject  unrelated  to  thermoacoustics, 
but,  nevertheless  of  interest  to  ONR,  namely  the  propagation  of  pressure  waves  through  layers  of  bubbly 
liquid,  both  with  the  purpose  of  attenuating  the  wave  and  to  generate  low-frequency  sound  by  efficient 
parametric  amplification.  Since  this  topic  is  outside  the  focus  of  interest  of  the  rest  of  this  final  report,  we 
will  only  give  a  concise  presentation.  Details  are  available  in  a  paper  submitted  to  J.  Acoust.  Soc.  Am.. 

A  schematic  representation  of  the  situation  considered  here  is  shown  in  Fig.  4.1:  a  one-dimensional  layer 
of  liquid  containing  gas  bubbles  is  located  between  x  —  0  and  x  —  L  and  is  excited  by  a  plane  wave  normally 
incident  from  the  left.  As  a  result  of  this  excitation,  a  reflected  wave  at  the  left  of  the  layer  and  a  transmitted 
wave  at  the  right  are  generated. 
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Figure  4.1:  Schematic  representation  of  the  one-dimensional  bubble  layer  excited  by  a  normally  incident 
plane  wave  from  the  left. 
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4.1  Mathematical  model 


The  mathematical  model  of  the  bubbly  liquid  consists  of  the  continuity  equation 


1  dP  du  _  d£ 

picf  dt  dx  dt 


in  which  pi  and  c*  are  the  density  and  sound  speed  of  the  pure  liquid  and  P,  u  the  average  mixture  pressure 
and  velocity,  and  of  the  momentum  equation 


du 

PlHi 


+  £-■ 


(4.2) 


This  model  is  essentially  that  of  Kogarko  (1964)  and  van  Wijngaarden  (1968)  except  that,  as  pointed  out 
by  Caflisch  et  al.  (1985),  the  convective  term  of  the  material  derivative  can  be  omitted  due  to  the  assumed 
smallness  of  the  gas  volume  fraction  0;  additional  considerations  on  this  point  are  given  in  Watanabe  & 
Prosperetti  (1994),  and  further  applications  of  this  and  similar  models  can  be  found  e.g.  in  Zabolotskaya 
(1977),  Kuznetsov  et  al.  (1978),  Gasenko  et  al.  (1979),  Nigmatulin  (1991),  Akhatov  et  al.  (1994),  Naugol- 
nykh  &  Ostrovsky  (1998),  Colonius  et  al.  (2000),  and  many  others.  Buoyancy  effects  are  neglected  in  (4.2) 
due  to  the  smallness  of  the  acoustic  time  scale  compared  with  the  time  evolution  of  the  bubble  layer.  The 
volume  fraction  is  given  by 

0(x,t)  =  ^nR\x,t)n  (4.3) 

where  R(x,  t)  is  the  instantaneous  radius  of  the  bubbles  contained  in  a  small  volume  centered  around  x  and 
n  is  the  bubble  number  density.  In  the  same  assumption  «  1  under  which  (4.1)  and  (4.2)  hold,  the  bubble 
number  density  n  can  be  taken  as  independent  of  time;  for  simplicity,  we  further  assume  it  to  be  spatially 
uniform  and  we  also  assume  that,  at  equilibrium,  all  the  bubbles  have  the  same  radius. 

In  spite  of  its  appearance,  the  previous  model  retains  a  strong  nonlinearity  in  the  manner  in  which  R  is 
calculated.  Again  on  the  basis  of  the  smallness  of  0,  for  this  purpose  we  use  the  Rayleigh-Plesset  equation 
of  bubble  dynamics  (see  e.g.  Plesset  &  Prosperetti  1977;  Prosperetti  1991;  Feng  &  Leal  1997) 


R 


d2R 
dt 2 


1  /  2ir  4/i  OR  \ 

ft  P  “  ”  ~R  ~  ~R~dt  ) 


(4.4) 


Here  p  is  the  bubble  internal  pressure  (approximated  by  the  gas  pressure,  the  small  vapor  contribution  being 
neglected),  a  the  surface  tension  coefficient,  and  p  the  liquid  viscosity.  For  an  isolated  bubble,  the  ambient 
pressure  P  appearing  in  (4.4)  is  to  be  identified  with  the  pressure  at  the  location  of  the  bubble  if  the  bubble 
were  absent.  In  a  dilute  mixture  the  bubbles  are  subject  to  the  averaged  field  and  P  should  be  taken  as 
the  average  pressure  appearing  in  the  momentum  equation  (4.2)  (see  e.g.  Caflisch  et  al.  1985;  Zhang  & 
Prosperetti  1994).  As  before,  in  (4.4)  we  omit  the  convective  term  of  the  material  derivatives  of  R. 

In  order  to  close  the  system  a  relationship  between  the  gas  pressure  p  and  bubble  radius  R  is  needed. 
This  point  has  been  treated  at  length  in  earlier  papers  (Prosperetti  1991;  Watanabe  &  Prosperetti  1994). 
Suffice  it  to  say  that  we  approximate  the  gas  pressure  inside  each  bubble  as  spatially  uniform,  which  leads 


to 


dp 

dt 


3 

r,  ,  dT 

dR] 

R 

(1-1)*  ^ 

n  JP9t. 

(4.5) 


where  T  is  the  local  gas  temperature  to  be  found  from 


7  p  (dT  8T\ 
7  -  1  T  (  dt  +V  dr  ) 
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Figure  4.2:  Steady  state  shape  of  the  transmitted  wave  according  to  the  complete  model  (solid  line),  compared 
with  poly  tropic  (dotted  line)  and  quasi-equilibrium  (dashed  line)  models  with  kp  =  1;  upper  panel  Pine /Poo  —  0.01, 
lower  panel  Pinc/Poa  ~  0.7;  w/27r  —  1.347  kHz,  f3o  —  0.133%,  Ro  —  0.121  mm. 


In  these  equations  7  and  k  =  k(T )  are  the  ratio  of  specific  heats  and  thermal  conductivity  of  the  gas  and  r 
is  the  radial  coordinate  measured  from  the  center  of  the  bubble.  As  shown  in  Kamath  et  al.  (1993).  at  the 
surface  of  the  bubble,  a  suitable  boundary  condition  for  (4.6)  is 

T(R,t)  =  (4.8) 

where  Tx,  is  the  undisturbed  liquid  temperature.  It  should  be  noted  that,  inside  the  bubble  centered  at  x, 
the  temperature  field  T  depends  on  r  as  well  as  t  and  hence,  in  principle,  the  set  of  equations  (4.5)  to  (4.7) 
must  be  solved  at  all  spatial  location  in  the  layer. 

The  model  described  in  this  section  is  supplemented  by  the  conditions  of  continuity  of  pressure  and 
velocity  across  the  planes  delimiting  the  bubbly  layer. 

In  addition  to  the  model  just  described,  in  this  work  we  consider  also  two  simpler  versions,  one  in  which 
the  internal  pressure  is  given  by  a  polytropic  relation 


and  one  in  which  the  same  relation  is  used,  but  the  left-hand  side  of  the  Rayleigh-Plesset  equation,  describing 
the  effect  of  the  liquid  inertia  on  the  radial  motion  of  the  bubbles,  is  set  to  zero.  We  refer  to  these  two 
simplified  models  as  to  the  poly  tropic  and  quasi- equilibrium  model,  respectively. 

The  model  equations  are  solved  numerically  by  a  total-variation-diminishing  method  similar  to  the  one 
described  in  Chapter  1. 

4.2  Results 

Figure  4.2  shows  the  transmitted  wave  at  steady  state  normalized  by  the  amplitude  of  the  incident  wave  as 
a  function  of  the  dimensionless  time  £cj/27t  for  two  incident  amplitudes,  Pine /Poo  —  0-01  (upper  panel)  and 
Pind Poo  —  0.7  (lower  panel).  The  solid  lines  are  the  results  of  the  complete  model  outlined  in  the  previous 
section.  The  dotted  and  dashed  lines  are  the  results  for  the  polytropic  and  quasi-equilibrium  models  with  kp 
=  1,  respectively.  When  the  incident  amplitude  is  small  (upper  panel)  the  polytropic  and  quasi-equilibrium 
models  give  essentially  the  same  result,  which  is  expected  in  this  case  in  which  the  incident  frequency  is  much 
smaller  than  the  bubble  resonance  frequency.  The  complete  model  predicts  a  somewhat  smaller  amplitude 
and  a  small  phase  shift,  both  due  to  the  inclusion  of  thermal  dissipation  in  the  bubble  motion,  but  behaves 
otherwise  very  similarly. 

For  the  larger-amplitude  excitation  (lower  panel  of  Fig.  4.2)  all  three  models  predict  shock  formation,  but 
the  differences  among  them  are  more  pronounced.  The  complete  model  (solid  line)  shows  a  transmitted  wave 
with  slight  oscillations  near  the  maximum,  while  these  oscillations  are  highly  exaggerated  by  the  polytropic 
model  (dotted  line).  The  same  qualitative  difference  is  encountered  when  the  two  models  are  applied  to 
shock  waves  in  bubbly  liquids  (Watanabe  &  Prosperetti  1994).  Since  these  oscillations  are  a  consequence  of 
bubble  inertia,  their  absence  in  the  quasi-equilibrium  model  (dashed  line)  is  not  surprising.  The  reason  why 
this  feature  is  encountered  at  this  higher  amplitude  but  not  at  the  lower  amplitude  of  the  upper  panel  of 
Fig.  4.2  is  that  the  formation  of  the  shock  introduces  a  much  shorter  characteristic  time  scale  in  the  wave, 
which  is  not  too  far  from  the  resonant  period  of  the  bubbles.  Due  to  the  absence  of  thermal  effects,  the 
polytropic  model  is  less  dissipative  than  the  complete  one,  and  the  shock  time  scale  is  accordingly  shorter: 
the  smaller  damping  and  the  shorter  time  scale  combine  to  cause  the  prominent  oscillations  of  this  result. 
If  these  oscillations  are  averaged  out  in  the  mind’s  eye,  one  sees  a  substantial  similarity  among  the  different 
profiles  which  is  due  to  the  fact  that  the  underlying,  relatively  slowly- varying,  wave  structure  is  only  slightly 
damped  in  all  models.  In  any  event,  it  may  be  noted  that  the  high-frequency  oscillations  of  the  complete 
model  are  so  strongly  damped  that  the  quasi-equilibrium  model  ends  up  being  a  better  approximation  to 
the  actual  behavior  than  the  polytropic  model. 

Figure  4.3  shows  the  components  of  the  transmitted  (upper  panel)  and  reflected  (lower  panel)  waves  at 
the  incident  frequency  for  Pinc/ P^  =  0.05  (circles)  and  Pinc/Poo  -  0.7  (squares)  according  to  the  complete 
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Figure  4.3:  Pressure  amplitudes  of  the  component  at  the  incident  wave  frequency  u;  for  the  transmitted  (upper 
panel)  and  reflected  (lower  panel)  waves  as  functions  of  the  incident  frequency  w  normalized  by  the  first  quasi- 
equilibrium  linear  eigenfrequency  of  the  layer,  =  1.347  kHz.  The  circles  and  squares  are  the  numerical  results 

of  the  complete  model  for  amplitudes  Pine! Poo  =  0.05  and  Pine /Poo  =  0.7  respectively.  The  solid  curves  are  the 
linear  results.  The  dotted  lines  connecting  squares  are  only  guides  to  the  eye.  The  bubble  radius  is  1.21  mm. 
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Figure  4.4:  Relative  power  of  the  low-frequency  component  in  the  transmitted  (circles)  and  reflected  (squares) 
waves  for  dual-frequency  excitation  as  a  function  of  the  amplitude  Pine /Poo  of  the  incident  wave  components 
for  u/1)  =  9uji  —  27rx3.03  kHz,  u)^  =  IOcji  =  27tx3.37  kHz  (Poo  is  the  undisturbed  ambient  pressure).  The 
lines  are  only  guides  to  the  eye.  The  solid  and  dotted  lines  are  for  the  complete  model  while  the  dot-dash 
and  dashed  lines  are  for  the  quasi-equilibrium  model  with  k,p  =  1.  The  bubble  radius  is  50  fim. 


model  as  functions  of  the  normalized  incident  frequency  oj/wu  with  o;1/27t  =  1.347  kHz  as  before;  the 
dotted  lines  connect  the  symbols  as  an  aid  to  the  eye.  The  solid  line  shows  analytic  results  obtained  from 
a  linearization  of  the  governing  equations  and  is  in  excellent  agreement  with  the  numerical  results  for  the 
lower-amplitude  case.  Remarkably,  when  normalized  by  the  incident  pressure  as  here,  the  higher-amplitude 
results  are  only  slightly  different,  with  somewhat  less  pronounced  maxima  and  minima  and  a  slight  shift  to 
lower  frequencies. 

Another  interesting  case  is  that  in  which  the  incident  wave  contains  two  frequencies: 

Pinc  =  A  jeos  uj^(x/c  —  t)  +  cos  oj^2\xjc  —  i)j  ,  (4.10) 

where  c  is  the  speed  of  soundf  in  the  liquid.  In  particular,  since  the  eigenfrequencies  of  the  normal  modes 
of  a  layer  such  as  the  one  considered  here  are  all  multiples  of  the  lowest  eigenfrequency,  one  can  take  ca(1) 
and  o/2>  to  correspond  to  eigenmodes  and,  in  addition,  their  difference  a/2)  —  a/1)  also  to  coincide  with  an 
eigenmode.  It  may  be  expected  that,  in  this  way,  there  would  be  a  very  efficient  parametric  effect.  In  the 
example  that  follows  we  take  ca(1)  and  a/2)  equal  to  the  9th  and  10th  eigenfrequencies,  so  that  their  difference 
coincides  with  the  fundamental. 

Figure  4.4  shows  the  relative  power  of  the  low-frequency  signal  (Pn/Pinc)2  of  the  transmitted  (circles) 
and  reflected  (squares)  waves  as  a  function  of  the  amplitude  of  the  incident  wave  components  Pinc/Poo  when 
cjH)  =  9^  and  a/2)  =  10wi;  the  solid  and  dotted  lines  are  for  the  complete  model  while  the  dot-dash  and 
dashed  lines  are  for  the  quasi-equilibrium  model.  Here  the  bubble  radius  is  Ro  =  50  /un,  with  a  linear 
resonance  frequency  u)o/2n  =  55.86  kHz  and  kp  =  1.01.  It  can  be  seen  here  that  the  efficiency  of  conversion 
is  about  1%,  which  is  a  very  high  value. 
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Figure  4.5:  Steady  state  shape  of  transmitted  wave  produced  by  the  complete  model  at  time  t  =  8  x  2iz/u\  for 
dual- frequency  excitation  with  —  9coi  ~  27rx3.03  kHz,  =  lOan  —  27rx3.37  kHz  and  amplitude  Pine/ Poo  — 

0.5. 


Figure  4.5  shows  the  steady  state  temporal  transmitted  waveform  at  time  t  =  8  x  2'k/uj\  for  a/1)  =  9cu1? 
-  lOcJi,  Pinc /Poo  =  0.5,  for  50  ^m-radius  bubbles.  The  transmitted  wave  has  a  sawtooth  appearance, 
but  the  peaks  are  not  very  sharp,  which  is  evidence  of  the  stronger  damping  affecting  the  higher  frequencies. 
The  relative  power  of  the  difference-frequency  harmonic  in  the  transmitted  wave  in  this  example  is  7.6  x  10~3. 


4.3  Conclusions 

The  results  summarized  here  correct  and  expand  the  scope  of  the  earlier  ones  of  Druzhinin  et  al.  (1996). 
In  comparison  with  that  work,  we  have  used  a  more  realistic  description  of  the  bubble  behavior  which 
includes  the  effects  of  radial  inertia,  with  the  associated  dispersion,  and  of  the  gas  thermal  behavior,  with 
the  attendant  energy  losses.  As  expected,  several  details  of  the  predictions  of  the  earlier  model  are  modified 
by  these  effects,  although  the  basic  character  remains.  Thus,  we  find  a  propensity  for  shock- wave  generation 
in  the  bubbly  layer  which  gives  rise  to  a  transmitted  wave  with  a  sawtooth  character. 

We  have  also  studied  the  possibility  of  enhanced  low-frequency  difference-wave  generation  through  the 
exploitation  of  the  resonances  of  the  bubbly  layer.  We  have  concluded  that  a  difference- wave  power  of  the 
order  of  1%  of  the  incident  power  can  be  generated  using  incident  wave  amplitudes  of  less  than  100  kPa. 
Thus,  this  technique  for  the  parametric  generation  of  low-frequency  waves  may  have  practical  value. 

An  important  aspect  of  the  phenomenon,  confirmed  by  this  study,  is  that  operation  near  the  resonance 
frequency  of  the  individual  bubbles  is  detrimental  to  the  energy  conversion  efficiency  due  to  the  strong 
dissipation  of  the  bubble  motion  in  this  frequency  range. 

In  order  to  avoid  the  practical  difficulties  connected  with  the  generation  and  control  of  suitable  bubbles, 
it  might  be  expedient  to  apply  the  principle  described  in  this  paper  to  other  systems.  For  example,  one 
would  expect  similar  phenomena  to  occur  in  porous  water-like  (or  rubber-like)  media  in  which  the  shear 
modulus  is  small  and  plays  the  role  of  gas  compressibility  in  bubbles;  some  experiments  of  this  type  are 
reported  in  Belyaeva  &  Timanin  (1991).  Like  bubbles,  pores  in  such  media  provide  very  strong  nonlinearity 
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(see,  e.g.  Naugolnykh  &  Ostrovsky  1998  section  1.4).  At  the  same  time,  it  is  much  easier  to  have  small  and 
almost  equal-size  pores,  the  system  is  more  stable,  and  losses  are  typically  smaller. 
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Chapter  1 

Introduction  and  Summary 


This  part  of  the  ONR-funded  project  focused  primarily  on  the  development  of  two-dimensional,  unsteady, 
physical  model  and  direct  numerical  simulation  scheme  for  the  study  of  simplified  thermoacoustic  refrigera¬ 
tors. 

The  model  is  based  on  the  so-called  zero-Mach-number  approximation  of  the  compressible  conservation 
equations.  The  advantages  of  this  asymptotic-based  approach  is  that  it  enables  us  to  overcome  the  large 
scale  disparity  that  characterizes  thermo-acoustic  devices.  Specifically,  the  approach  combines  a  simplified 
representation  of  the  action  of  acoustic  waves,  with  detailed  multi-dimensional  resolution  of  the  flow  and 
temperature  fields  in  the  neighborhood  of  thermoacoustic  stacks  and  heat  exchangers. 

Based  on  this  attractive  physical  modeling  approach,  a  vorticity-based  numerical  scheme  was  developed 
for  efficient  simulation  of  thermal  and  mechanical  behavior  of  thermo- acoustic  refrigerators.  Specifically,  the 
numerical  scheme  combines  domain  decomposition  /  boundary  Green’s  function  techniques  with  FFT-based 
elliptic  equations  solvers,  which  results  in  highly  efficient  approach  for  the  simulation  of  thermo-acoustic 
devices. 

At  the  start  of  this  project,  we  started  with  a  reduced  formulation  based  on  an  adiabatic  flow  assumption 
and  a  physical  model  incorporating  a  single  stack  with  no  heat  exchangers.  This  initial  formulation  was 
subsequently  relaxed  and  the  model  underwent  several  levels  of  refinement.  Various  milestones  have  been 
summarized  in  a  number  of  archival  publications  [1,  2,  3,  4,  5,  6,  7],  and  the  refinements  culminated  in 
Dr.  Besnoin’s  Ph.D.  thesis  [8]. 

In  Chapter  2,  a  summary  of  the  most  refined  version  of  the  physical  model  is  presented.  Chapter  3  then 
outlines  the  corresponding  numerical  scheme.  We  conclude  in  Chapter  4  with  a  summary  of  major  findings. 


Chapter  2 


Physical  Formulation 


The  physical  formulation  has  been  specifically  adapted  to  the  simplified  thermoacoustic  refrigerator  geometry 
illustrated  in  figure  2.1.  The  governing  equations  are  derived  by  following  the  procedure  described  in  [6] 
which  is  summarized  below.  First,  the  Navier-Stokes  equations  are  made  dimensionless  using  an  appropriate 
choice  of  normalizing  parameters.  The  gas  dynamic  variables  are  then  expanded  in  terms  of  the  Mach  number 
and  the  corresponding  asymptotic  series  is  substituted  for  each  variable,  and  the  leading  order  equations  are 
derived. 


Acoustic  Driver  Resonance  Tube  Rigid  End  of 


Figure  2.1:  Schematic  illustration  of  the  thermoacoustic  refrigerator  (top)  with  a  magnified  view  of  the 
computational  domain  (bottom). 
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2.1  Governing  equations 

The  formulation  is  based  on  the  assumptions  that:  (a)  the  amplitude  of  the  pressure  oscillations  in  the 
inner  region  is  small  compared  to  the  background  pressure;  (b)  the  flow  field  is  two-dimensional;  (c)  the 
stack  and  heat  exchangers  are  non-intrusive,  meaning  that  they  don’t  affect  the  flow  in  the  “outer  region”, 
which  remains  quasi-ID:  (d)  the  stack  and  heat  exchangers  have  dimensions  that  are  much  smaller  than  the 
wavelength  of  the  acoustic  standing  wave,  so  that  the  inner  computational  domain  can  be  assumed  to  be 
acoustically  compact  [2,  6];  and  (e)  the  gas  is  Newtonian  and  has  constant  viscosity  and  thermal  conductivity. 

Under  these  assumptions,  the  Mach  number  (M  =  uref/c)  is  small  and  a  low  Mach  number  expansion 
of  each  gas  dynamic  in  an  asymptotic  series  in  terms  of  e  =  7 M2  is  appropriate.  Here,  7  is  the  specific  heat 
ratio,  uref  is  the  reference  velocity  and  cref  is  the  sonic  velocity.  By  expanding  the  equation  variables  in  an 
asymptotic  series  in  terms  of  M,  and  by  suitably  combining  the  state,  continuity  and  energy  equations,  the 
normalized  Navier-Stokes  equations  governing  the  motion  of  the  gas  are  reduced  to  [9,  2]: 

1.  The  equation  of  state: 


P  =  pT 


(2.1) 


2.  The  equation  for  the  evolution  of  the  pressure: 


VdP 
7  dt 
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which  can  also  be  written  as: 
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3.  The  Helmholtz  decomposition  is  used  to  decompose  the  velocity  field  into  an  irrotational  part  V</>  and 
a  divergence  free  rotational  part  V  x  (ipk): 


u 

V2ip 

V2<p 


V<j>  +  V  x  (ipk) 
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4.  The  vorticity  transport  equation: 


<K 

dt 
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5.  The  energy  equation  for  the  fluid: 
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dT 

dt 


—pu  -  VT  + 
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Pe  Re 


(2.4) 

(2.5) 

(2.6) 


(2.7) 


(2.8) 


Here  u  =  (u,  v )  =  u/uref  is  the  velocity  vector,  ( the  vorticity,  <f>  the  velocity  potential,  ip  the  streamfunc- 
tion,  k  is  the  right-handed  normal  to  the  plane  of  motion,  P  the  thermodynamic  pressure,  T  the  temperature, 
p  the  density,  $  the  viscous  dissipation  function  [10],  D/Dt  =  d/dt  +  u  ■  V  the  material  derivative,  7  the 
specific  heat  ratio,  while  D  denotes  the  computational  domain,  dD  its  boundary,  and  n  the  outer  normal 
at  dD.  Re  =  ureflref  /u  is  the  Reynolds  number,  v  is  the  kinematic  viscosity.  Pe  =  prefCpureflref/k  is  the 
Peclet  number  for  the  gas,  k  is  the  thermal  conductivity  of  the  gas  and  C'p  is  the  specific  heat  at  constant 
pressure  of  the  gas.  Ec  =  pref /(CpTref)  is  the  Eckert  number.  fref,  pref ,  lref,  are  the  reference  tempera¬ 
ture,  fluid  density  and  length,  respectively. 
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The  temperature  distribution  within  the  plates  is  governed  by  the  unsteady  heat  conduction  equation: 

\  Xi'2T  (2  9) 

dt~PesVTs ’  1  } 

where  Pes  =  Q.H2 /ds)  is  the  Peclet  number  for  the  solid.  Finally,  the  hot  and  cold  heat  exchangers  are 
treated  as  isothermal.  Their  respective  temperatures  7),  and  Tc  are  imposed. 

One  of  the  key  features  of  the  present  model  is  that  it  enables  us  to  isolate  a  small  neighborhood  of  the 
stack  and  heat  exchangers.  Within  this  region,  which  is  assumed  to  be  “acoustically  compact”,  the  model 
ignores  elastic  effects  but  retains  bulk  compressibility  and  convective  nonlinearities  in  a  stratified  medium. 
This  enables  us  to  avoid  sonic  CFL  limitations  and  thus  efficiently  analyze  the  details  of  the  2D  flow. 

2.2  Normalization 

The  variables  are  normalized  with  respect  to  the  appropriate  combination  of  angular  frequency  fi  of  the 
acoustic  standing  wave,  stack  plate  spatial  periodicity  H,  mean  pressure  Pm,  mean  temperature  fm,  and 
mean  density  of  the  fluid  pm.  Thus,  we  have: 

•  ^ref  — 

®  Ire/  —  ^ 

•  Pref  =  Pm 

•  Tref  =  Tm 

•  Pref  “  Pm 

The  choice  of  H  and  Q.  -1  as  length  and  time  scales  leads  to  the  reference  velocity  scale  ure/  =  DF. 

2.2.1  Flow,  fluid  and  material  characteristic  parameters 

The  normalization  process  leads  to  the  definition  of  the  following  dimensionless  parameters: 

•  Reynolds  number,  Re  =  flH2/D 

•  Eckert  number,  Ec  =  pmj (CpTm) 

•  Gas  Peclet  number  Pe  =  f \H2  /a 

•  Solid  Peclet  number,  Pes  =  QH2/as,  which  describes  the  thermal  properties  of  the  plate  material. 

Here  as  and  a  =  k/(pCp)  denote  the  thermal  diffusivities  of  the  solid  and  the  gas,  respectively.  In  addition  to 
these  four  dimensionless  numbers  which  appear  during  the  derivation  of  the  governing  equation,  the  following 
parameters  are  relevant  in  the  study  of  the  thermally  stratified  flows  in  thermoacoustic  refrigerators: 

•  The  viscous  penetration  depth  5V  =  y2f>/G.  Within  a  distance  5V  from  the  solid  surfaces,  viscous 
shear  forces  cause  gradients  in  oscillating  velocities. 

•  The  thermal  penetration  depth  Sk  =  y25//fi.  The  thermal  penetration  depth  is  a  scale  for  the 
distance  that  heat  diffuses  through  the  fluid  in  an  acoustic  period.  In  conventional  thermoacoustic 
engines,  8k  ranges  from  0  (0.1  mm)  to  0  (1  mm),  and  is  orders  or  magnitude  smaller  than  the  acoustic 
wavelength  A.  Gas  particles  located  further  than  8k  from  the  solid  surfaces  undergo  essentially  adiabatic 
oscillations  and  do  not  contribute  to  the  thermoacoustic  effects. 

•  The  particle  displacement  amplitude  Rp  =  ua/0,  where  ua  =  cDrsin(kx)/j,  is  the  local  velocity  am¬ 
plitude  derived  using  linear  acoustics  in  the  resonance  tube.  In  the  literature,  the  particle  displacement 
is  sometimes  referred  to  as  the  tidal  displacement,  which  is  equivalent  to  2 Rp. 
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•  The  Prandtl  number  Pr  =  p.Cp/k  =  (/)„/Sk)2.  The  Prandtl  number  is  a  dimensionless  measure  of  the 
ratio  of  viscous  to  thermal  effects.  In  the  case  of  monatomic  ideal  gases,  Pr  ~  2/3. 

•  The  acoustic  Reynolds  number  Re„  =  2ua/V$&.  The  acoustic  Reynolds  number  is  a  dynamic  Reynolds 
number  based  on  the  acoustic  wave  velocity  ua-  Note  that  Re„  —  i  Re  R2.  Rea  can  be  used  to  determine 
whether  the  Stokes  boundary  layers  are  stable  or  undergo  transition  to  turbulence  [11]. 

2.2.2  Geometrical  parameters 

The  geometry  of  the  resonance  tube,  stack  and  heat  exchangers  (Fig.  2.2)  is  specified  in  terms  of: 

•  E,  the  stack  plate  spatial  periodicity,  or  centerline  separation  distance. 

•  d,  the  thickness  of  plate  and  heat  exchangers.  In  the  idealized  refrigerator  depicted  in  figure  2.2,  the 
plates  and  heat  exchangers  are  assumed  to  have  the  same  thickness. 

•  h  is  the  spacing  between  two  consecutive  plates  or  heat  exchangers.  We  have  H  =  d+h. 

•  The  blockage  ratio,  BR  =  h/H,  which  is  used  to  study  the  blockage  effect  of  the  stacks  on  the  flow. 

•  Lh  and  Lc,  which  denote  the  length  of  the  hot  and  cold  heat  exchangers,  respectively.  When  Lh  and 
Lc  have  the  same  value,  both  parameters  are  synthetizised  using  Lh/C. 

•  g,  the  gap  width  that  separates  the  thermoacoustic  stack  from  the  stacks  of  heat  exchangers.  The 
distance  between  the  stack  plates  and  the  cold  heat  exchangers  and  the  distance  between  the  stack 
plates  and  the  hot  heat  exchangers  can  be  different.  However,  in  all  the  computations,  both  distances 
are  considered  equal,  which  justifies  the  absence  of  subscripts  in  the  notation. 

•  Dhe-dD  is  the  distance  which  determines  the  location  of  the  acoustic  matching  surfaces. 

•  kx  is  the  dimensionless  wave-number  based  on  the  non-dimensional  wave-number  k  and  the  distance 
of  the  stack  from  the  rigid  end  of  the  tube  x. 

All  the  parameters  are  listed  in  the  nomenclature. 


2.3  Boundary  Conditions 

The  simulations  are  based  on  the  assumption  that  the  thermoacoustic  stack  is  composed  of  a  large  number 
of  identical,  equi-spaced  plates  and  that  the  hot  and  cold  heat  exchangers  have  the  same  number  of  plates 
and  the  same  plate  spacing  H  as  the  stack  (fig.  2.2).  We  model  this  setting  by  restricting  the  simulations  to 
a  single  period  and  imposing  periodic  boundary  conditions  in  the  direction  normal  to  the  acoustic  axis. 
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Figure  2".2:  Schematic  illustration  of  the  computational  domain,  showing  the  dimensions  of  the  stack  and 
heat  exchanger  plates,  the  gap  g  between  the  stack  and  heat  exchangers,  the  channel  height  H  and  the 
location  of  the  matching  boundaries. 
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As  shown  in  figure  2.3,  the  boundary  of  the  computational  domain  is  divided  into  8  sets,  dD  ~  (jjLi  T*. 
In  the  cross-stream  direction,  the  computational  domain  is  limited  by  the  sets  r2,  F3,  r5  and  T6,  where 
periodic  boundary  conditions  axe  imposed.  The  inlet /outlet  surfaces  Ti  and  T4  limit  the  computational 
domain  in  the  “streamwise”  direction.  On  these  surfaces,  the  solution  matches  the  wave  representation  in 
the  “outer  region”,  which  is  assumed  to  be  dominated  by  a  quasi-ID,  idealized  acoustic  standing  wave.  In 
the  following,  we  refer  to  these  computational  surfaces  as  matching  boundaries. 


Hot  HX  Cold  HX 


Figure  2.3:  Schematic  illustration  of  the  computational  domain  and  the  limiting  boundaries. 


2.3.1  Condition  at  Matching  Surfaces 

The  action  of  the  standing  wave  is  represented  by  prescribing  oscillations  of  the  thermodynamic  pressure, 
P,  and  the  mean  flow  velocity.  Here,  the  integral  constraint  in  (Eq.  2.3)  is  used  to  compute  the  oscillating 
velocities  at  the  ends  of  the  computational  domain,  which  in  turn  are  used  to  impose  the  boundary  conditions 
for  the  vorticity,  streamfunction  and  velocity  potential  on  Fi  and  IV  The  computed  inlet  and  outlet  velocities 
?jj  and  U2  are  obtained  by  adequately  combining  equations  (2.10)  and  (2.11). 

The  difference  between  the  velocities  is  computed  using  Eq.  (2.3),  the  left  hand  side  of  which  can  be 
re-written  as: 

u-2  -ui  =  —  u  ■  n dA  (2.10) 

A  JdD 

The  second  constraint  is  obtained  assuming  that  the  average  of  the  inlet  and  outlet  flow  velocities  matches 
the  prescribed  instantaneous  mean  flow  velocity  ua  at  the  location  of  the  stack.  ua  is  obtained  by  using 
linear  acoustics  in  the  resonance  tube.  This  leads  to:  - 


^(u2  +  «i)  =  ua  =  ^  sin(fcr)  cos(f) 


(2.11) 


12]: 


Thermal  boundary  conditions  at  the  matching  surfaces  correspond  to  the  vanishing  flux  conditions  [3,  4, 


dT 

—  =0  on  Ti  and  T4 
oy 


(2.12) 


These  conditions  are  motivated  by  the  fact  that  as  one  moves  away  from  the  stack  and  heat  exchangers  the 
thermoacoustic  effect  is  expected  to  essentially  vanish,  so  that  the  axial  thermal  gradients  should  decrease 
as  well. 

2.3.2  Boundary  conditions  in  the  cross-stream  direction 

In  the  fluid  and  the  stack  plate,  the  periodicity  boundary  conditions  on  surfaces  T2  to  T5  are  implemented 
by  imposing: 

4>(y,t)  lr2  =  0(2/>*)lr5  (2-13) 

ip(y,t)  lr2  =  'i/’(y-*)lr5  (2-14) 

T(y,t)  |ra  =  T(y,t)  |r#  (2-15) 
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and 


Ts(y,t)\ri  =  Ts(y,t)\Vg  (2.16) 

As  noted  in  [2],  the  imposition  of  periodicity  conditions  on  the  potential  (2.13)  and  streamfunction  (2.14)  is 
desirable,  since  it  prevents  the  development  of  a  mean  flow  along  the  y  direction,  and  thus  avoids  potential 
difficulties  associated  with  imposing  periodicity  conditions  on  the  velocity  itself. 

r9  is  an  internal  boundary,  which  represents  the  surface  of  the  plate;  IV  and  Ts  represent  the  surfaces 
of  the  cold  and  hot  heat  exchangers  respectively.  At  these  3  surfaces,  no  slip  and  zero-normal  velocity 
conditions  are  used  and  are  implemented  by  imposing, 


d(j> 

=  0 

(2.17) 

dn 

=  0 

(2.18) 

d'ljj 

<90 

(2.19) 

dn 

dp 

where  n  is  the  normal  to  the  plate  surface  and  p  —  n  x  k  is  the  unit  vector  tangent  to  the  plate  surface.  Note 
that,  as  discussed  in  2.3.1,  the  particular  choice  of  velocity  potential  and  streamfunction  boundary  conditions 
corresponds  to  letting  the  potential  velocity  carry  the  entire  volume  flux  through  the  domain.  Thus,  the 
streamfunction  distribution  does  not  induce  a  mean  flow  and  describes  a  non-linear  vortical  perturbation  to 
an  otherwise  ideal,  irrotational  flow. 

2.3.3  Thermal  boundary  conditions  at  the  solid  surfaces 

The  heat  exchangers  are  treated  as  isothermal  solid  bodies  with  constant  hot  temperature  Th,  and  cold 
temperature  Tc,  which  are  imposed  on  Ts  and  Ty  respectively. 

Meanwhile,  on  T9,  the  conduction  solution  within  the  stack  plates  is  coupled  to  the  energy  equation  for 
the  gas.  We  require  that  the  temperature  solutions  and  the  corresponding  conduction  fluxes  within  the  gas 
and  solid  coincide  at  the  plate  surface.  Thus,  we  impose: 

T(q,t)  -  Ts(q,t)  (2.20) 

£<•••>  -  1 :«•*>  (2-2I> 

where  q  is  a  coordinate  variable  along  T9  and  n  =  ks/k  is  the  ratio  of  the  thermal  conductivity  in  the  solid 
to  the  thermal  conductivity  of  the  gas. 
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Chapter  3 

Numerical  Scheme 


The  numerical  simulation  of  the  model  equations  is  based  on  a  finite  difference  methodology.  Spatial  deriva¬ 
tives  are  approximated  using  second-order  centered  differences,  and  time  integration  is  based  on  the  third- 
order  Adams-Bashforth  scheme. 

The  construction  of  the  scheme  can  be  decomposed  into  two  parts.  First,  the  model  variables  are  updated, 
following  a  procedure  which  varies  according  to  the  type  of  representation  for  the  acoustic  standing  wave 
discussed  in  2.3.1.  Both  procedures  are  developed  in  section  3.2.  The  second  step  consists  in  reconstructing 
the  velocity  fields  from  the  vorticity  and  divergence  fields.  The  reconstruction  of  the  velocity  involves  the 
solution  of  two  elliptic  equations  which  are  inverted  using  fast  Poisson  solvers.  The  fast  solvers  are  based 
on  domain-decomposition  /  boundary  Green’s  functions  techniques,  the  construction  of  which  is  discussed 
in  3.3. 


3.1  Discretization  and  notation 

The  numerical  implementation  developed  in  [6]  was  extended  in  order  to  accommodate  the  presence  of  heat 
exchangers.  The  discretization  also  carries  out  a  wider  variety  of  boundary  conditions  at  the  matching 
surfaces.  The  model  equations  are  discretized  using  a  finite  difference  method.  Specifically,  the  model 
variables  are  discretized  on  a  rectangular  mesh  with  size  Ax  and  A y  in  the  x  and  y  directions  (Fig.  3.1). 
In  the  cross-stream  direction,  NX  is  used  to  denote  the  number  of  cells  between  P2  and  T5  in  domains  £>i, 
Ds,  D5  and  D7,  whereas  NX2  denotes  the  number  of  cells  between  two  consecutive  plate  or  heat  exchanger 
surfaces.  In  the  streamwise  direction,  NY,  indicates  the  number  of  cells  in  the  corresponding  domain  DL. 


HotHX  -P.+.T  ■u.v.tv  Cold  HX 


Figure  3.1:  Schematic  illustration  of  the  computational  domain  mesh  grid  and  of  the  internal  boundaries 

Velocity,  streamfunction  and  vorticity  distributions  are  discretized  on  the  edges  of  the  computational 
cells.  Density,  temperature  and  velocity  potential  are  discretized  on  the  cell  centers  as  shown  on  Fig  3.1. 
Time  integration  of  the  system  of  discrete  evolution  equations  is  performed  using  the  explicit,  third-order 
Adams-Bashforth  scheme. 

Spatial  derivatives  are  discretized  using  second-order  differences.  On  a  rectangular  grid,  a  standard 
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center-difference  approximation  to  the  gradient  and  Laplacian  operators  respectively  yields: 

Vfc/&  = 

and 


fn 

tn  J  *+*  J 


fn  fn  _  fn 

I  i— 1J  ,  fij+l  JiJ— 1 
2Ax  2A  y 


fi+u  ~  2/&  +  fi-ij 
Ax2 


V2  fn  _  't-t-id  -  *-,j  •  *—±*j  i 

hi  id  “  h 


fn  __  O  f  n  f  n 

J id+ 1  _  lid-l 

Ay'2 


(3-1) 


(3.2) 


The  material  derivative  at  time  nAt  is  approximated  using: 


Dt 


du™  ■ 

-Al  +  uljVhulj 

3<i-K71  +  ^2 

2At 


(3-3) 


i.e.,  the  temporal  derivative  is  estimated  using  a  second-order  backward  difference  while  the  convective  term 
is  computed  using  a  standard  centered-difference  formula. 

As  shown  in  figure  3.1,  the  computational  domain  for  the  flow  field  can  be  divided  into  seven  rectangular 
internal  sub-domains  £>i  to  D7,  bounded  in  the  streamwise  direction  by  the  matching  boundaries  Ti  and  r4, 
and  by  six  internal  boundaries  to  r(y/).  Boundaries  F(/)  to  r(v/)  are  computational  boundaries  used 
in  the  numerical  solution  of  elliptic  equations  (sections  3.3.1  and  3.3.2). 


3.2  Time  stepping 

Implementation  of  the  numerical  scheme  to  the  simulation  of  the  unsteady  2D  flow  is  based  on  successive 
repetition  of  the  following  steps. 

1.  Update  the  pressure,  the  pressure  derivative  and  the  mean  flow  velocity: 

Pn+1  =  1  -f  Dr  cos (kx)  sin[(n  +  1)A t] 

(jip  \  n~\~l 

—  j  =  Dr  cos  (kx)  cos[(n  +  l)At] 


Dr 

7  M 


sin  (kx)  cos[(n  H-  l)t] 


(3-4) 

(3.5) 

(3.6) 


Compute  the  right  hand  side  (rhs3)  of  equation  2.3  : 


+ 


Ec 

Re 


L*"‘dv] 


(3.7) 

J  u 

Combine  equations  2.10  and  2.11  so  as  to  update  the  velocities  at  the  matching  surfaces  rd  and  F>: 

(3.8) 


<+1  =  <+1  -  \rhsff 


<+1  =<+t  +  Irhs^1 


(3-9) 

'  ' 

2.  Integrate  equation  (2.7)  in  order  to  determine  the  vorticity  values  in  the  interior  of  the  domain  at  the 
new  time  step.  We  set: 


,n—h 


2  vy  n—k  n  n- 

C#1  =<&  +  **!>[  -  V/,  x  («”Jfc  x  C?Jk)  X  ~Dt 

k=0 


rid 


(3.10) 
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3.  Integrate  the  energy  equation  (2.8)  in  order  to  determine  the  new  temperature  distribution  at  interior 
points  in  the  fluid  domain.  We  use: 


rpn+ 1  _  rnn  ,  Ay. 


-  «57*-vfc25y*  + 


7  —  1  dP 


<n—k 


k—0 


ipfjk 


dt 


+ 


Pep?” 


-r-r'2  rpn—k 
:n-k  Vh1i,j 


+ 


—$n~k' 
Re  «  J 


(3-11) 


4.  Integrate  the  heat  conduction  equation  (2.9)  in  order  to  update  the  temperature  field  at  interior  points 
within  the  solid.  Thus,  we  set: 


rp  in-fl 

lsHJ 


.„  At 

-  +  - - 

Pes 


I>  117 


k- 0 


(3.12) 


5.  Determine  the  temperature  distribution  at  solid  surfaces  by  requiring  that  the  heat  flux  be  continuous 
there,  i.e.  by  enforcing  (eq.  2.21): 


l  surf 


in+ 1  _ 


(4Tj 


in+1 

h-ij 


. + 4«^++P)  -  cr.iak + ^ij) 

3f  1  +  k) 


(3.13) 


6.  Determine  the  new  density  field  in  the  fluid  from  the  equation  of  state  (2.1).  We  set: 

pn+l 


Pn+1  = 


(3.14) 


7.  Determine  the  new  streamfunction  distribution  and  velocity  potential  by  inverting  the  discrete  Poisson 
equations  corresponding  to  equations  (2.5)  and  (2.6).  This  step  is  discussed  in  sections  3.3.1  and  3.3.2. 

8.  Reconstruct  the  velocity  field  from  the  divergence  and  vorticity  fields  (eq.  2.4)  using  the  centered- 
difference  approximations  to  the  gradient  operator  described  in  3.1. 

9.  And  finally,  update  the  vorticity  at  the  surface  of  the  heat  exchangers  and  plate  by  imposing  no-slip 
conditions  at  the  solid  surfaces.  For  more  details,  see  [4]. 

Note  that  steps  1-6,  8  and  9  only  involve  straightforward  algebraic  manipulation.  On  the  other  hand, 
step  7  requires  the  inversion  of  two  discrete  elliptic  operators,  which  accounts  for  most  of  the  CPU  time  in 
the  computations.  The  inversion  of  the  streamfunction  Poisson  equation  and  determination  of  the  potential 
distribution  are  performed  using  the  fast  domain-decomposition  solver  discussed  in  sections  3.3.1  and  3.3.2. 
The  efficiency  of  both  solvers  was  examined  in  [2].  Comments  on  the  advantages  of  the  inversion  technique 
can  be  found  in  [6]. 


3.3  Domain  decomposition 

The  inversion  procedures  for  the  determination  of  the  streamfunction  distribution  and  the  velocity  potential 
are  summarized  below.  The  procedures  were  first  developed  by  Worlikar  and  Knio  in  [2,  6,  4]  and  are  adapted 
to  accommodate  the  inversion  of  the  Laplacian  operator  in  the  computational  domain  shown  in  figure  3.1. 
Here,  the  technique  is  extended  to  the  decomposition  of  the  computational  domain  into  7  sub-domains. 

3.3.1  Streamfunction  Distribution 

The  streamfunction  distribution  at  the  new  time  level  is  governed  by  the  following  Poisson  equation, 


with  homogeneous  Dirichlet  boundary  conditions  on  Fi,  F4  and  Tg,  and  periodicity  conditions  on  r2/r5. 
The  inversion  is  performed  using  a  domain  decomposition  technique  which  subdivides  the  computational 
domain  into  seven  non-overlapping  regions,  D  =  (j[=1  £>;•  As  sketched  in  figures  2.3  and  3.1,  Di  is  bounded 
by  the  surfaces  r i .  T?,  Fg,  Ts  and  F -F .  Do  by  F-F ,  Fg  and  Do,  by  F.,,  Fg  and  r9, 

Di  by  r<m>,  r9  and  r</y>,  D5  by  r</v\  r^r,,  r5,  rT  and  r9j  D6  by  r(y),  rr  and  and  Dr  by 

r(v/>,r2,  r5,  r7  and  r4. 

The  streamfunction  distribution  within  individual  sub-domains  is  decoupled  using  a  boundary  Green’s 
function  technique  which  is  based  on  constructing,  in  a  pre-processing  step,  M  elementary  solutions  in  each 
sub-domain,  in  Du  l  =  1,2, 3, 4, 5, 6, 7,  m  =  1,...,M,  where  M  is  the  number  of  internal  boundary 
points  on  rG\r(r/),r(//J),F(/vr),  r(V)  and  F(v/>,  {zm}* f=1.  The  boundary  Green’s  functions  are  the 
solutions  of  the  following  elementary  problems: 


V3A=  0  inF»' 


(3.16) 


l  =  1, ..., 7,  j  =  1,...,M,  with  homogeneous  Dirichlet  boundary  conditions  on  r4,  r4,  r7,  r8  and  F9, 
periodicity  conditions  on  r2  and  r5,  and  the  following  boundary  conditions  on  the  internal  boundaries 


=  /  1 
i  1  0 


if  x  =  Zj 
otherwise 


(3.17) 


The  individual  streamfunction  solutions  at  the  new  time  level  are  expressed  as  a  linear  combination  of  the 
boundary  Green’s  and  of  auxiliary  solutions.  We  thus  write: 


ibLn+l  =  yl'n+l  + 


M 

E 

m— 1 


’!m  m 


1  =  1,...,  7, 


(3.18) 


where  the  auxiliary  solutions  xj/'n+l  are  determined  by  inverting: 

V|$',n+1  =  -^'n+1  (3.19) 

with  homogeneous  Dirichlet  conditions  on  F'  and  T",  periodicity  conditions  on  F2  and  r5,  and  homogeneous 
Dirichlet  conditions  on  Fi,  F4  and  T7.  The  unknown  coefficients  7j"+1  are  determined  by  requiring  that  the 
normal  derivative  of  the  streamfunction  at  internal  boundaries  -approximated  within  individual  sub-domains 
using  second-order  one-sided  derivatives-  be  continuous  at  -  r(v/F  This  results  in  an  M  x  M  system 
of  linear  equations,  whose  corresponding  influence  matrix  is  denoted  by  A.  Its  inverse,  A-1,  is  determined 
by  Gaussian  elimination  in  a  pre-processing  step,  and  then  fed  to  the  computations.  The  advantage  of  the 
domain  decomposition  /  boundary  Green’s  function  approach  is  that  it  decouples  the  inversion  procedures 
in  elementary  sub-domains  and  facilitates  the  implementation  of  fast  solvers.  In  the  computations,  we  take 
advantage  of  the  fact  that  the  sub-domains  Dx  -  D7  are  rectangular,  and  of  the  nature  of  the  boundary 
conditions  on  horizontal  surfaces.  Specifically,  in  Dt,  D3,  D5  and  D7  we  exploit  of  the  periodicity  conditions 
on  r2/F5  by  using  a  discrete  FFT  of  the  governing  equations  [13,  14].  In  Do,  D\  and  D§,  we  take  advantage 
of  the  homogeneous  Dirichlet  conditions  on  the  streamfunction  by  implementing  a  fast  sine  transform.  These 
transformations  lead  to  decoupled  tridiagonal  systems  for  the  Fourier/ sine  coefficients,  which  are  efficiently 
determined  using  the  Thomas  algorithm  [13,  14].  Thus,  the  streamfunction  distribution  is  determined  using 
a  procedure  that  requires  essentially  0 (N  log  N)  operations. 


3.3.2  Potential  Distribution 


The  velocity  potential  is  governed  by  the  following  Poisson  equation: 


hi 


1  dPn+1 
7  dt. 


+  Pe  ,l  fR»  y 


Re 


(3.20) 


with  periodicity  conditions  on  T2  and  r5,  non-homogeneous  Neumann  conditions  on  r4  and  T4  given  by 
equations  (2.10,2.11),  and  a  homogeneous  Neumann  condition  on  the  solid  surfaces,  F7,  Fg  and  Fg.  Fast 
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solution  of  the  above  system  is  performed  using  a  domain  decomposition  technique  which  was  developed  by 
Worlikar  and  Knio  [4]  and  which  is  summarized  below. 

The  construction  of  the  velocity  potential  solver  is  more  involved  than  that  of  the  streamfunction  solver, 
in  large  part  due  to  the  Neumann  boundary  conditions  imposed  on  the  solid  surface,  Tj,  Tg  and  T9  and 
the  acoustic  matching  surfaces,  Tx  and  F4.  One  difficulty  is  that  the  use  of  Dirichlet  conditions  on  the 
coupling  surfaces  T(/)  -  would  lead  to  a  mixed  Dirichlet/Neumann  problem  on  T^/Ts  -  F^^/Tr, 
which  would  greatly  complicate  (or  prevent)  implementation  of  a  fast  Fourier  solver.  An  additional  concern 
is  that  the  potential  distribution  is  only  determined  up  to  an  additive,  arbitrary  constant;  thus,  the  coupling 
of  solutions  across  sub-domain  boundaries  must  carefully  account  for  this  feature.  Also,  sufficient  care  must 
be  exercised  so  that  the  compatibility  condition  for  the  original  problem  (3.20),  which  is  guaranteed  in  light 
of  (2.2),  remains  satisfied.  In  order  to  address  these  issues,  we  first  note  that  in  the  implementation  of  the 
domain  decomposition  technique  it  is  only  necessary  to  determine,  within  individual  sub-domains,  potential 
distributions  that  (i)  satisfy  compatibility  conditions  locally  over  each  sub-domain,  and  (ii)  yield  velocity 
predictions  that  are  continuous  across  computational  surfaces  joining  different  sub-domains  (F(/l  to  F(v/)). 
This  is  the  case  because  it  is  only  the  gradient  of  the  velocity  potential  that  enters  in  the  formulation  of  the 
physical  model.  Based  on  this  observation,  we  proceed  as  follows.  We  use  the  same  domain  sub-division 
used  in  section  3.3.1  and  express  the  solution  within  individual  sub-domains  as: 

M’ 

<j)i'n+i  =  ^n+i + 1  =  i>->7>  (3*21) 

k=  1 


where  M'  is  the  total  number  of  segment  centers  lying  on  T W  -  T^VI\  {zrk}^f=l.  Note  that  M’  =  M  -F  6, 
and  that  the  collocation  points  along  the  interface  differ  from  those  used  for  the  streamfunction,  i.e.  zk  ^  zk 
for  all  k.  For  clarity  of  the  presentation,  we  assume  that  the  collocation  points  zl  are  indexed  consecutively 
such  that  zk  £  for  1  <  k  <  AT ,  zfk  £  F^7*  for  Ny  +  1  <  k  <  2 TV7,  zlk  £  for  2 Ny  4- 1  <  A;  <  3AT7, 

zlk  £  T(/V)  for  3 JV7  + 1  <  k  <  4AT7,  z'k  £  r(v’>  for  4 JV7  +  1  <  k  <  5AT7  and  zfk  £  T{VI)  for  5AT7  + 1  <  k  <  M'. 
The  elementary  boundary  Green’s  functions  are  elementary  solutions  of  the  following  non-homogeneous 
Neumann  problems: 

i=e',  1  =  1,...,  7,  (3.22) 


where 


(3.23) 


and  \Di\  is  the  volume  of  domain  Dj.  In  conjunction  with  (3.22),  we  use  periodicity  conditions  on  F>  and 
r5,  homogeneous  Neumann  conditions  on  rx,  T4,  r7,  Ts  and  r9,  and  the  following  Neumann  conditions. 


d$‘ 


dy 

d*i 


dy 


-  =  1  if  x  =  z‘k 
—  0  otherwise 


(3.24) 


on  rW  -  T^VIl  Note  that  since  the  $lk  are  solutions  of  Neumann  problems  they  are  only  determined  up 
to  an  additive  arbitrary  constant,  and  that  the  compatibility  condition  for  (3.22)  is  automatically  satisfied 
following  the  definition  of  el  (3.23).  The  potential  distributions  <?’n+\  l  =  1, ...,  7,  are  the  numerical  solutions 
of: 

V|0z,n+1  =  e*’n+1  (3.25) 

with  periodicity  conditions  on  T2  and  T5,  homogeneous  Neumann  conditions  on  r7,  F8  and  T9,  and  the 
following  inhomogeneous  Neumann  conditions  on  Ti,  T4,  -  T^VI\ 
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(3.31) 

(3.32) 

(3.33) 


respectively.  Note  that  the  solvability  conditions  for  (3.25)  in  D1  -  D6  are  automatically  satisfied  by  con¬ 
struction  of  the  non-homogeneous  Neumann  conditions  (3.28)  to  (3.33),  respectively,  while  the  compatibility 
condition  for  (3.25)  in  D7  follows  from  the  compatibility  condition  of  the  original  problem  (3.20).  Also  note 
that,  by  construction,  the  combination  of  boundary  Green’s  functions  and  auxiliary  solutions  (3.21)  auto¬ 
matically  satisfies  the  continuity  of  the  normal  derivative  of  0  at  r(/)  -  r(v/h  Thus,  the  obvious  next  step 
would  be  to  determine  the  unknown  coefficients  cr”+1  by  requiring  continuity  of  the  potential  itself  at  r(/) 
-  T^.  However,  such  an  exercise  would  be  complicated  since  both  the  auxiliary  solutions  0/,n+1  and  the 
boundary  Green’s  functions  are  determined  up  to  an  arbitrary  additive  constant.  To  avoid  this  difficulty, 
while  at  the  same  time  accounting  for  the  indeterminacy  of  0,  we  first  observe  that  if  a  reference  value  for 
the  potential  in  each  subdomain  is  correctly  selected  at  one  point  along  the  interface,  then  continuity  of  the 
potential  at  the  remaining  points  would  follow  from  continuity  of  the  tangential  derivative.  The  selection  of 
the  reference  point  is  of  course  arbitrary;  for  convenience,  we  choose  the  first  collocation  points  on  r(/)  - 
r<v/)  ^  reference  points,  and  so  require  the  continuity  of  the  tangential  derivative  at  the  remaining  points. 


This  results  in: 
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i.e.  a  linear  system  of  M'  -  6  equations  in  the  M'  unknown  Green’s  function  coefficients.  To  complete  the  sys¬ 
tem,  two  additional  equations  are  used  which  require  that  the  overall  source  term  (divergence)  corresponding 
to  (3.21)  coincides  to  that  of  the  original  (undecomposed)  problem  (3.20).  Thus,  we  impose: 

Nt 

^<r”+1=0  (3.35) 

fc= i 


53  tf+1  =  o  (3.36) 

k=Nr+l 

The  combined  influence  matrix,  B,  corresponding  to  (3.34),  (3.35)  and  (3.36)  is  inverted  in  a  pre-processing 
step,  and  its  inverse  is  stored  for  use  in  the  unsteady  simulations. 


3.4  Concluding  remarks 

/ 

With  the  present  procedure,  the  unknown  coefficients  of  the  boundary  Green’s  functions  are  uniquely  de¬ 
termined.  The  end  result  is  a  potential  distribution  that  satisfies  the  original  governing  equation  within 
individual  subdomains,  and  guarantees  continuity  of  the  normal  and  tangential  derivatives  of  the  potential 
along  computational  boundaries  separating  the  subdomains.  Since  the  physical  model  only  involves  the 
gradient  of  the  velocity  potential,  continuity  of  the  potential  itself  does  not  need  to  be  guaranteed.  The  pri¬ 
mary  advantage  of  the  present  approach  is  that  it  facilitates  the  implementation  of  fast  solvers  for  individual 
Neumann  problems  in  (3.22)  and  (3.25).  In  the  present  computations,  we  exploit  the  special  structure  of  the 
problems  by  using  a  fast  Fourier  solver  Di,  D%,  D-t ,  and  D-  and  a  fast  cosine  transform  in  IDo .  and  Dq . 
Consequently,  the  overall  algorithm  is  also  0(N  log  N).  Note  that  although  the  streamfunction  and  velocity 
potential  are  discretized  using  different  grids,  the  Fourier-based  techniques  used  in  the  fast  inversion  of  the 
corresponding  Poisson  problems  are  compatible,  i.e.  they  rely  on  the  same  number  of  points.  The  advantage 
of  using  a  staggered  grid  for  <t>  is  that  the  potentially-delicate  task  of  imposing  Neumann  conditions  at  the 
plate  corners  is  avoided.  In  the  simulations,  computational  grids  are  used  that  have  a  relatively  large  number 
of  points.  As  a  consequence,  it  is  impractical  to  store  the  boundary  Green’s  functions  for  the  streamfunction 
and  the  potential.  To  overcome  this  problem,  the  inversion  of  the  Poisson  equations  for  the  streamfunction 
and  the  velocity  potential  is  performed  twice.  The  first  inversion  enables  us  to  determine  the  correct  condi¬ 
tions  at  the  computational  boundaries  T(/)  -  r(V  /),  while  the  second  inversion  yields  the  complete  solution. 
Despite  the  additional  effort,  the  utilization  of  the  fast  solver  remains  quite  attractive  [6]. 
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Chapter  4 

Summary  of  Major  Results 


Direct  numerical  simulations  have  been  extensively  applied  in  a  fundamental  investigation  of  thermoacoustic 
refrigerators,  with  specific  attention  to  thermal  and  mechanical  behavior  of  the  stack  and  heat  exchangers. 
These  analyses  fell  within  the  following  four  areas: 

•  Quantitative  visualization  of  the  flow  field  in  the  neighborhood  of  thermoacoustic  stacks  [1,  2,  6]; 

•  Analysis  of  mechanical  energy  losses  [4]; 

•  Validation  studies  and  fundamental  analysis  of  thermal  stack  behavior  [3,  8];  and, 

•  Optimization  of  heat  exchanger  configuration  [7,  8]. 

Highlights  from  these  studies  are  briefly  discussed  below. 


4.1  Flow  Visualization 

As  briefly  outlined  above,  a  detailed  study  of  the  unsteady  flow  field  was  carried  out  in  [1,  2].  For  brevity, 
for  one  configuation  from  the  case  study  in  [2]  is  described  here. 

The  details  of  the  flow-field  are  illustrated  for  a  stack  with  a  blockage  ratio  BR  =  0.666  and  a  plate 
length  parameter  L/H  =  1.32.  The  Reynolds  number  and  the  drive  ratio  are  Re  =  2132  and  Dr  =  1.0%, 
respectively.  The  stack  is  kept  at  a  position  mid-way  between  the  velocity  node  and  the  anti-node  (kx  = 
37r/4). 

The  model  is  allowed  to  run  till  a  quasi-steady  periodic  state  is  reached,  i.e.  the  flow  shows  the  same 
variation  over  consecutive  acoustic  cycles.  Figure  4.1  shows  the  streamlines  plotted  for  one  acoustic  cycle. 
The  frames  are  generated  one-eighth  of  a  cycle  apart  and  the  arrows  at  the  left  of  the  frame  indicate  the 
magnitude  and  direction  of  the  mean  velocity  at  the  center  of  the  stack.  The  first  frame  is  plotted  at  the 
start  of  the  acoustic  cycle  when  the  acoustic  velocity  at  the  center  is  zero.  The  flow  starts  accelerating 
upwards  till  the  velocity  is  maximum  in  frame  (c).  The  flow  then  starts  decelerating  till  the  velocity  vanishes 
in  frame  (e).  The  fluid  then  starts  accelerating  in  the  downward  direction,  reaches  a  maximum  downward 
velocity  in  frame  (g)  and  starts  decelerating  again.  This  sequence  is  repeated  after  frame  (h). 

As  the  flow  starts  accelerating  upward  (frame  (b)),  recirculation  bubbles  are  formed  just  at  the  entrance 
of  the  channel.  At  the  same  time,  an  eddy  is  shed  from  the  exit  of  the  channel.  An  eddy  at  the  entrance 
which  was  formed  in  the  earlier  cycle,  gets  crushed  against  the  forward  face  of  the  channel  and  slides  over 
the  recirculation  eddy  formed  at  the  entrance  of  the  channel  and  migrates  towards  the  centerline  forming  a 
vortex  bubble  near  the  centerline  of  the  channel  which  then  travels  with  the  flow.  As  the  flow  accelerates 
upward,  the  recirculation  zone  at  the  entrance  of  the  channel  moves  upward  within  the  channel.  Meanwhile, 
the  separation  bubbles  at  the  exit  of  the  channel  continues  to  intensify  as  does  the  separation  bubble  in 
the  channel  near  the  entrance.  As  the  flow  is  reversed  (frames  e-h),  these  processes  are  reversed  -  i.e.  a 
recirculation  eddy  is  formed  behind  the  bottom  end  of  the  channel,  whereas  a  recirculation  zone  is  formed 
at  the  throat  at  the  top  end  of  the  channel. 
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Figure  4.1:  Evolution  of  streamfunction  distribution  for  a  stack  with  h/H  =  0.66,  L/H  =  1.32,  Dr  =  1.0%. 
The  frames  are  generated  at  times  (a)  t  =  327t/4,  (b)t  =  33tt/4,  (c)t  =  347r/4,  (d)  t  =  35tt/4,  (e)  t  —  367r/4, 
(f)  t  =  377r/4,  (g)  t  =  387r/4  and  (h)t  =  397r/4. 


As  illustrated  in  figure  4.1,  our  analysis  reveals  that  the  flow  field  is  strongly  modulated  by  the  motion 
of  complex  vortical  structures  which  form  due  to  the  shedding  of  Stokes  layers,  and  by  the  impingement  of 
these  structures  on  solid  surfaces.  As  discussed  below,  these  motions  substantially  affect  the  heat  transfer 
rates  and  thermal  performance  of  the  refrigerator. 


4.2  Mechanical  Energy  Losses 

Extensive  computations  based  on  a  simplified  version  of  the  low-Mach-number  number  were  carried  out  in  [5] 
in  order  to  analyze  the  fluid  force  balance  and  mechanical  energy  losses  in  a  thermoacoustic  stack  at  low  to 
moderate  drive  ratios. 

In  particular,  this  study  has  tested  the  hypothesis  that  the  net  hydrodynamical  effect  on  the  momentum 
balance  of  a  stack  of  narrowly  spaced  plates  in  an  otherwise  essentially  one-dimensional  resonant  acoustic 
field  can  be  represented  by  effective  impedance  laws  similar  to  those  used  in  linear  theory.  To  this  end,  the 
analysis  assumed  a  purely  sinusoidal  input  of  mean  mass  flux  and  weak  (acoustic)  compressibility  to  drive  an 
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inelastic  two-dimensional,  weakly  compressible  flow  around  a  stack  of  rectangular  plates.  The  forces  imposed 
by  the  plates  on  the  surrounding  flow  were  consequently  analyzed. 

The  force  balance  was  first  used  to  show  that  the  response  of  the  stack  on  the  outer  sinusoidal  driving  is 
sinusoidal  as  well,  up  to  errors  of  no  more  than  2%  in  all  the  cases  considered.  Explicit  simplified  formulae 
of  the  plate  forces  were  then  derived  based  on  the  linear  theory.  The  simplest  version  assumes  Poiseuille 
flow  between  the  plates  at  all  times,  a  more  sophisticated  one  accounts  for  the  dependence  of  the  velocity 
profile  on  inertia  and  viscous  effects.  By  comparison  with  the  force  balance  from  the  detailed  numerical 
simulation  good  agreement  was  obtained  for  the  elaborate  theory,  but  unsatisfactory  agreement  was  found 
for  the  simplified  Poiseuille-flow  model. 

From  the  force  estimates  explicit  impedance  laws  were  further  obtained  that  can  be  used  in  quasi- 
one-dimensional  simulations  as  approximate  models  of  the  effect  of  the  thermoacoustic  stack  on  the  outer 
long  wave  acoustics.  Comparison  of  the  semi-heuristic  estimates  for  the  impedance  with  the  results  of  the 
computations  showed  that  for  plate  length  to  plate  spacing  ratios  larger  than  about  L/H  =  8  and  for 
Reynolds  numbers  in  the  range  of  20  <  Re  <  550,  the  linear  theory  yields  reliable  estimates  of  the  effective 
impedance  of  the  stack.  A  distinctive  feature  of  the  impedance  calculations  in  [5]  is  that  they  account  for 
finite  thickness  effects.  This  is  a  crucial  aspect  since  the  2D  flow  simulations  indicate  that  the  pressure  drag 
exerted  by  the  plate  has  a  comparable  magnitude  to  the  friction  force,  even  for  stacks  having  long  plates 
and  blockage  ratios  as  large  as  0.9. 

4.3  Validation  and  Analysis 

Detailed  computations  have  also  been  performed  in  order  to  establish  the  validity  of  the  model  predictions  [6, 
3,  8],  and  to  analyze  the  limitations  of  linearized  predictions. 

A  sample  of  these  validation  studies  are  provided  in  figure  4.2,  which  contrasts  computed  predictions  and 
and  experimental  measurements  of  the  steady  state  temperature  difference  across  an  unloaded  thermoacoustic 
stack.  Good  agreement  is  evident  between  the  experimental  results  and  the  computations  at  all  the  values 
of  the  drive  ratio  considered.  The  computations  are  in  excellent  agreement  with  the  predictions  of  the 
theory  at  low  drive  ratio,  Dr  <  1%.  However,  as  Dr  increases,  the  predictions  of  the  linear  theory  deviate 
significantly  from  both  the  experimental  measurements  and  the  computations.  At  moderate  drive  ratios, 
the  deviation  can  be  substantial,  and  exceeds  25%  at  a  drive  ratio  of  0.02.  In  addition,  as  previously  noted 
in  [15],  the  theoretical  estimates  tend  to  overpredict  AT  at  all  stack  positions,  except  in  the  neighborhood  of 
the  velocity  node  where  the  thermoacoustic  effect  naturally  vanishes.  Note  that  the  computed  stack  position 
where  AT  peaks,  kx  ~  2.8,  coincides  with  the  corresponding  experimental  measurement.  On  the  other  hand, 
the  simulations  and  experiments  are  in  slight  disagreement  with  the  linear  theory,  which  predicts  that  the 
peak  temperature  difference  occurs  at  kx  ~  2.7. 

The  distribution  of  the  energy  flux  density  in  the  neighborhood  of  the  cold  end  of  the  stack  is  plotted  in 
figure  4.3  for  two  values  of  the  drive  ratio.  At  Dr  =  0.5%,  figure  4.3  shows  that  as  they  leave  the  cold  end  of 
the  stack  the  energy  pathlines  rapidly  curve  around  the  corner  and  then  follow  an  essentially  straight  path 
towards  the  hot  end.  Meanwhile,  at  Dr  =  2%,  figure  4.3  shows  that  as  they  exit  the  end  and  side  of  the 
plate  the  energy  pathlines  first  contour  a  large  well-defined  recirculation  region  before  they  enter  into  the 
channel.  Thus,  for  the  presently  considered  conditions,  the  mean  energy  path  between  the  ends  of  the  stack 
significantly  increases  with  increasing  drive  ratio.  This  phenomenon,  and  to  some  extent  the  occurrence  of 
a  relatively  important  2D  motion  within  the  gap,  are  the  major  effects  of  the  increase  in  Dr  and  appear  to 
be  at  the  origin  of  the  departure  between  the  theory  and  the  computations. 
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Figure  4.2:  Top:  Variation  of  temperature  difference  across  the  stack  with  drive  ratio.  The  stack  is  located 
at  kx  =  37t/4.  Bottom:  Variation  of  temperature  difference  across  the  stack  with  position  at  Dr  =  2%. 
Experimental  data  are  adapted  from  [15],  and  correspond  to  TAC  #  3.  The  stack  plates  are  a  0.1905-mm 
thick,  laminate  of  stainless  steel  and  fiberglass,  with  an  effective  thermal  conductivity  ks  =  5.76  W/m-°K. 
The  stack  is  6.85-mm  long  with  a  half  gap  yo  =  h/2  =  0.669  mm.  Thus,  the  stack  has  a  blockage  ratio 
h/H  =  0.875  and  a  plate  length  parameter  of  L/H  —  4.48.  The  physical  dimensions  of  the  stack  are 
approximated  in  the  computations  by  using  a  blockage  ratio  h/H  =  0.837  and  a  plate  length  parameter 
L/H  —  4.5.  Also  shown  are  the  predictions  of  the  linear  theory. 
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Computations  have  also  been  performed  in  [8]  to  validate  velocity  field  predictions.  A  sample  of  this 
validation  study  is  provided  in  figures  4.4  and  4.5,  which  contrast  computations  and  measurements  of  the 
instantaneous  velocity  field  in  the  neighborhood  of  a  stack  edge. 

As  summarized  in  Table  4.1,  two  stack  configurations  are  considered  in  the  experiments  [16].  In  both 
cases,  the  center  of  the  stack  is  located  midway  between  the  velocity  node  and  the  velocity  anti-node,  i.e.  at 
a  non-dimensional  wave-number  kx  =  3n/4.  In  addition  to  the  geometrical  parameters  of  the  stack,  the 
computations  use  as  input  the  thermal  properties  of  the  stack  plates  and  of  air,  as  well  as  the  mean  pressure 
and  temperature,  respectively  Pm  =  105  N/m2  and  Tm  =  300°K.  The  experimentally-observed  frequency  of 
the  resonant  standing  wave,  /  =  200  Hz  is  also  used  as  input.  Meanwhile,  the  amplitude  of  the  standing 
wave  is  expressed  in  terms  of  the  drive  ratio,  Dr ,  which  expresses  the  ratio  of  the  acoustic  pressure  amplitude 
to  the  mean  pressure  (Table  4.1). 

The  two  stack  configurations  are  geometrically  similar,  and  are  nearly  equal  in  length.  However,  the 
thickness  and  spacing  of  the  stack  plates  are  significantly  smaller  in  configuration  A  than  they  are  in  con¬ 
figuration  B.  These  differences  are  also  highlighted  by  providing  the  ratios  h/5v  and  d/5v,  where  8V  is  the 
viscous  penetration  depth  [17].  Thus,  d  is  substantially  larger  than  5V  in  configuration  A,  but  is  close  to  5V  in 
configuration  B.  Accordingly,  we  shall  refer  to  these  two  stacks  as  thick-plate  and  thin-plate  configurations, 
respectively. 


d  (mm) 

h  (mm) 

Lp  (mm) 

Dr  (%) 

h/8v 

d/  deltav 

Configuration  A 

1.0 

2.0 

25.8 

0.5 

13. 

6.5 

Configuration  B 

0.15 

1.0 

24.0 

1.5 

6.7 

1.0 

Table  4.1:  Geometric  and  flow  parameters  for  configurations  A  and  B. 


Figure  4.4  depicts  two  instantaneous  realizations  of  the  velocity  and  vorticity  fields  for  configuration 
A.  Shown  on  the  top  row  are  experimental  results  obtained  at  to  +  T/4  and  to  +  3T/8;  the  bottom  row 
depicts  computational  predictions  obtained  at  tp  +  ST/16  and  tp  +  5T/16.  Thus,  in  both  the  experiments 
and  computations  the  two  frames  are  T/8  apart.  Note  that  the  frame  of  the  computational  results  has  also 
been  restrict  to  match  the  4.2  x  4.2  mm2  window  of  the  PIV  measurements  [16]. 

For  the  present  configuration,  both  the  experimental  and  computational  results  reveal  the  presence  of 
concentrated  vortices  near  the  edge  of  the  plate.  Also  evident  is  the  presence  of  both  signs  of  vorticity,  both 
inside  the  channel  and  in  the  open  region.  The  generation  of  alternating  layers  of  vorticity  within  the  channel 
is  not  surprising,  since  the  ratio  of  channel  height  to  viscous  penetration  depth  is  large  and  the  development 
of  Stokes  layers  within  the  channel  is  accordingly  expected.  As  can  be  observed  in  Figure  4.4,  the  vorticity 
distribution  outside  the  channel  can  also  exhibit  complex  structure,  especially  when  separated  vortices  are 
driven  back  into  the  channel  and  impinge  on  the  edge  of  the  plates. 

Comparison  of  the  experimental  and  computational  results  in  Figure  4.4  reveals  a  close  correspondence 
between  the  predictions.  In  particular,  at  both  phases,  very  good  agreement  can  be  observed  between  peak 
vorticity  values  as  well  as  the  size  of  the  recirculating  regions. 

Figure  4.5  provides  computed  and  experimental  results  obtained  for  configuration  B.  Shown  on  the  left 
column  are  computed  velocity  and  vorticity  distribution  at  three  selected  phases  within  the  cycle;  plotted 
on  the  right  columns  are  contours  of  the  experimental  vorticity  distribution  at  the  corresponding  times. 
The  figure  is  generated  in  a  similar  fashion  as  in  Fig.  4.4,  i.e.  the  computational  test  section  is  restricted  so 
that  it  matches  the  experimental  window.  In  the  present  case,  the  PIV  measurements  are  performed  in  a 
2.9  x  2.3  mm2  window,  also  located  near  the  cold  end  of  the  stack. 

The  present  predictions  are  thus  in  sharp  contrast  with  earlier  observations.  Specifically,  in  the  present 
case  the  results  do  not  show  the  formation  of  well-defined  eddies.  The  vorticity  distribution  in  the  wake 
exhibits  the  presence  of  elongated  layers  that  extend  well  outside  the  channel.  However,  significant  rollup  of 
these  layers  is  not  observed,  although  evidence  of  very  weak  recirculating  motion  near  the  plate  corners  can 
be  detected  in  the  last  computational  frame. 

Note  that  the  top  and  bottom  frames  in  Fig.  4.5  correspond  to  tp  and  tp+T/2,  times  at  which  the  acoustic 
pressure  is  maximal  and  minimal,  respectively.  At  these  phases,  the  streamwise  extent  of  the  vortical  wake,  is 
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found  to  be  maximal  and  the  direction  of  the  flow  in  the  middle  of  the  channel  is  opposite  to  that  of  the  flow 
along  the  plate  surface.  As  noted  by  Duffourd  [16],  these  observation  are  in  good  qualitative  agreement  with 
the  velocity  profiles  derived  using  the  linear  theory  [17,  18]  for  the  same  configuration.  The  computed  and 
measured  velocities,  which  are  in  good  agreement  with  each  other,  are  slightly  larger  than  those  predicted 
by  the  theory.  As  discussed  in  [2],  this  discrepancy  is  most  likely  due  to  blockage  effects  which  axe  ignored 
by  the  linear  theory  [2]. 


Figure  4.4:  Instantaneous  velocity  (vectors)  and  vorticity  (contours)  fields  around  the  cold  end  of  the  stack 
in  configuration  A.  Top  row:  experimental  PIV  measurements;  bottom  row:  computations.  The  times  at 
which  the  frames  are  generated  are  indicated  in  the  labels. 
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4.4  Optimization  of  Heat  Exchangers 

A  detailed  study  has  been  recently  conducted  in  [8]  of  the  effects  of  geometrical  and  operating  parameters 
on  the  performance  of  the  device.  Specifically,  a  parametric  study  has  been  performed  of  the  effect  on  the 
cooling  load  Qm,f->c  of  drive  ratio  (2%  <  Dr  <  8%),  temperature  difference  between  the  hot  and  cold 
heat  exchangers  (6  I<  <  AT  <21  K),  gap  width  g  between  the  heat  exchangers  and  the  stack  plates,  heat 
exchanger  length  Lc,  position  of  the  stack  in  the  resonance  tube  kx.  thickness  of  the  plates  d,  and  plate 
spacing  h. 

Analysis  of  the  computations  revealed  that  the  heat  exchanger  length,  gap  width,  drive  ratio,  temperature 
difference  between  the  heat  exchangers,  position  of  the  stack  in  the  resonance  tube,  heat  exchanger  and  plate 
spacing  have  a  combined  effect  on  the  performance  of  the  device.  More  specifically,  the  cooling  load  has 
been  found  to  peak  at  a  well-defined  combination  of  the  gap  width  and  heat  exchanger  length. 

Two  parameters,  namely  the  optimal  heat  exchanger  length  and  the  optimal  gap  width,  were  indepen¬ 
dently  determined.  The  optimal  length  varied  only  with  the  gap  width  and  with  the  plate  separation  distance; 
it  did  not  depend  on  any  of  the  other  parameters  studied.  On  the  other  hand,  the  optimal  gap  width  was 
independent  of  the  drive  ratio  only.  It  became  wider  when  (a)  the  temperature  difference  between  the  heat 
exchangers  was  increased,  (b)  the  plates  were  stacked  closer  together.  It  also  varied  with  the  stack  location 
in  the  tube. 

The  power  density  peaked  when  the  plates  were  separated  by  a  distance  about  2.5  times  the  thermal 
penetration  depth,  and  when  the  thickness  of  the  plates  wass  slightly  smaller  than  the  thermal  penetration 
depth.  With  thinner  or  thicker  plates,  the  power  density  dropped  sharply. 

Detailed  analysis  has  also  been  performed  in  [8]  of  the  impact  of  proper  optimization  of  heat  exchanger 
parameters.  In  particular,  it  has  been  shown  that  the  cooling  load  fluctuates  above  90%  of  its  maximal 
value  whenever  the  length  of  the  heat  exchangers  fell  within  0.8  <  Lh/J(2.RV  +  S/.)  <  2.  For  shorter  heat 
exchangers,  the  cooling  dropped  sharply. 

In  addition,  the  optimization  study  also  showed  that  the  gap  width  between  the  heat  exchangers  the 
stack  plates  should  be  accurately  selected.  The  introduction  of  a  finite  gap  between  the  stack  plates  and 
heat  exchangers  can  lead  to  a  significant  increase  in  the  cooling  load.  Specifically,  for  a  large  temperature 
difference  between  the  heat  exchangers,  optimization  of  the  gap  width  was  shown  to  have  a  dramatic  impact 
on  the  cooling  load  and  the  thermal  performance  of  the  device. 
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Appendix 


Nomenclature 

Roman 

BR 
Dr 

D/Dt 
Ec 
H 
La 
Lp 

Lg 

Lhjc 

P 
Pe 
Pe5 

Pr 

Qc—±f 
R 
Re 
Rea 

RP 
T 

V 
d 

9 
h 

kx 

P 
n 
t 

tp 

u 

Greek 

<3>  Viscous  dissipation 

a  Thermal  diffusivity  k/(pCp) 

dD  Boundary  of  the  computational  domain 

(j)  '■  Velocity  potential 

7  Specific  heat  ratio 

k  Ratio  of  gas  thermal  conductivity  to  plate  thermal  conductivity 


Blockage  ratio  (=  (H/(H  —  d)) 

Drive  ratio  (=  Po/Pm) 

Material  derivative  (=  d/dt  -f  u  •  V) 

Eckert  number  (~  Pm/(Cpfm)) 

Stack  centerline  spacing 
Normalizing  parameter  (=  2 Rp  +  5k) 

Length  of  the  stack  plates 

Length  of  the  gap  between  the  stack  and  each  heat  exchanger 
Length  of  the  hot/cold  heat  exchangers 
Thermodynamic  pressure 
Peclet  number  (=  pmQH2Cp/k) 

Peclet  number  for  the  plate  (=  Q.H2/as) 

Prandtl  number  (=  jlCp/k) 

Heat  transfer  from  the  cold  heat  exchanger  to  the  gas 
Specific  gas  constant 

Reynolds  number  (=  &H2/i>)  _ 

Acoustic  amplitude  Reynolds  number  (=  ua/V&&) 

Particle  displacement  parameter  (=  ua/QH) 

Temperature 

Volume  of  the  computational  domain 
Stack  and  heat  exchanger  thickness 
Spacing  between  stack  and  heat  exchangers 
Stack  gap  (=  2aj0) 

Dimensionless  wave  number  (=  flx/c) 

Hydrodynamic  pressure 

Outward  normal  to  the  boundary  of  the  computational  domain 
Time 

Time  at  which  the  pressure  rate  of  change  is  minimal 
Velocity  vector 
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n)  Streamfunction 

p  Density 

£  Vorticity 


Subscripts 

a  Acoustic 

c  Cold  heat  exchanger 

/  Fluid 

h  Hot  heat  exchanger 

m  Mean 

s  Solid 


Superscripts 

Dimensional  quantity 


Dimensional  variables 


cp 

k 

Po 

c 

n 

sv 

4 


V 

A 

Ur  <2 


-1  0-1  K”1) 


Isobaric  specific  heat  (J  kg”1  K"1) 

Thermal  conductivity  (J  m”1  s“ 

Pressure  amplitude  (N/m2) 

Speed  of  sound  (m/s) 

Angular  frequency  of  the  standing  wave  (rad/s) 

Viscous  penetration  depth  in  the  fluid  (m)  (=  yj 2v/Q) 

Thermal  penetration  depth  in  the  fluid  (m)  (=  y2o//0) 

Nominal  plate  thicloiess  (m) 

Kinematic  viscosity  (N-m/s) 

Wavelength  of  the  standing  wave  (m) 
local  velocity  amplitude  (=  cDrsin(kx)/ 7) 


Mathematic  Operators  and  Indexes 

rp 

V  Nabla  Operator,  =  ^ ,  ■§;) 

o  Dyadik  Product 

Scalar  Product 
A  t  Time  step 

Ax  Mesh  size  in  the  x-Direction 

A y  Mesh  size  in  the  y-Direction 
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1.  SUMMARY 

Long-term  research  objective:  The  long-term  objective  of  this  component  of  the  study  is  the 
experimental  investigation  of  heat  transfer  in  the  thermoacoustic  stack  and  heat  exchangers  with 
the  purpose  of  developing  better  models  and  designs  of  heat  exchangers. 

Technical  and  research  objectives:  The  heat  exchangers  of  a  thermoacoustic  refrigerator 
deliver  the  heat  absorbed  in  the  refrigerated  volume  to  the  cold  side  of  the  thermoacoustic  stack 
and  absorb  the  heat  transported  from  the  cold  to  the  hot  side  of  the  stack,  which  is  then  released 
to  the  ambient.  Improving  their  performance  is  key  to  improving  the  overall  performance  of 
thermoacoustic  refrigerators.  Only  few  studies  reported  in  the  literature  focus  on  the  modeling 
and  design  of  heat  exchangers  for  thermoacoustic  devices. 

Technical  approach:  The  approach  taken  on  our  research  is  to  first  visualize  and  measure  the 
high-speed  oscillating  temperature  fields  in  the  thermo  acoustic  stack  using  holographic 
interferometry  combined  with  high-speed  cinematography.  Results  of  the  heat  transfer 
measurements,  specifically  the  measured  heat  transfer  coefficients  were  used  as  input  data  for  the 
models  that  were  developed  to  analyze  the  performance  of  thermoacoustic  heat  exchangers. 

Publications:  During  the  past  years  6  papers  were  published  in  peer-reviewed  journals  (J1-J6, 
attached),  one  paper  (J7)  was  accepted  and  two  are  in  preparation  (J8  and  J9).  One  invited  book 
chapter  on  heat  transfer  in  oscillatory  flow  was  completed  (B2,  attached)  as  well  as  one  on  heat 
transfer  enhancement  of  heat  exchangers  (Bl).  A  contribution  was  made  to  the  Dictionary  of 
Acoustics  (B3).  Six  papers  were  published  in  proceedings  of  conferences  (C1-C6)  as  well  as  nine 
abstracts  in  conference  proceedings  (C7-C15).  The  research  has  resulted  in  one  dissertation  (Dl), 
two  Master’s  degrees  and  7  European  diploma  degrees  equivalent  to  a  Master’s  degree  in 
collaboration  with  universities  from  Europe. 
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Presentations:  The  work  has  been  presented  at  14  conferences  and  at  31  invited  lectures  in  the 
USA,  Switzerland,  Turkey,  France  and  Germany.  Two  invited  lectures  will  be  given  at  the 
NATO  Advanced  Study  Institute  on  Low-temperature  and  Cryogenic  Refrigeration  to  be  held 
this  summer  in  Cesme,  Turkey.  One  of  the  objectives  of  the  research  program  was  to  familiarize 
the  engineering  and  especially  heat  transfer  community  with  thermoacoustic  refrigeration.  This 
goal  was  accomplished  by  presenting  the  results  of  this  research  at  mechanical  engineering 
conferences  and  by  giving  seminars  at  research  labs,  institutions  and  universities  in  the  United 
States  and  Europe. 

Organization  of  conference  sessions:  A  session  on  Emerging  Refrigeration  Technologies  and 
one  on  Transport  Phenomena  in  Oscillatory  Flows  was  organized  by  Cila  Herman  and  Omar 
Knio  at  the  National  Heat  Transfer  Conference  in  Baltimore,  August  10-12 ,  1997.  Seven 
papers  were  submitted  for  the  sessions  and  four  of  the  presentations  focused  on  thermoacoustics. 
The  papers  were  published  in  the  proceedings  of  the  conference,  AIChE  Symposium  Series, 
Volume  93, 1997,  Mohamed  S.  El-Genk,  Volume  Editor.  Two  session  with  a  similar  topics  were 
organized  at  the  National  Heat  Transfer  Conference  in  Albuquerue,  New  Mexico ,  Aug .  15-19, 
1999 .  Cila  Herman  and  Martin  Wetzel  organized  the  Special  Session  on  Thermoacoustics  (2 
sessions,  6  invited  lectures,  16  papers)  at  the  137th  Meeting  of  the  Acoustical  Society  of  America 
and  the  2nd  Convention  of  the  European  Acoustics  Association  held  March  13-17,  1999  in 
Berlin,  Germany. 

Impact/Navy  relevance:  The  Navy  is  interested  in  efficient,  environmentally  safe  large-scale 
and  small-scale  refrigeration  systems,  for  example  for  air-conditioning,  cooling  of  food  or 
electronic  equipment,  etc.  The  development  of  more  efficient  heat  exchangers  can  lead  to 
smaller,  lighter  and  more  efficient  thermoacoustic  refrigerators,  which  is  of  particular  importance 
in  systems  where  size  and  weight  is  critical.  Systematic  design  optimization  tools  are  essential  in 
commercial  applications.  Miniature  thermoacoustic  refrigerators  (MEMS)  are  a  potential 
solution  for  a  range  of  practical  applications. 

Technology  transfer:  The  two  design  optimization  algorithms  developed  for  the  project  were 
used  to  evaluate  the  performance  of  the  miniature  thermoacoustic  refrigerator  for  electronic 
cooling  applications  and  improve  its  design  for  a  project  carried  out  with  Rockwell  Science 
Corporation,  CA.  A  series  of  lectures  was  given  at  the  Universite  Pierre  et  Marie  Curie  (Paris  6), 
Ecole  Centrale  Lyon  and  Ecole  Centrale  Paris  in  France  as  well  as  the  University  of  Hannover  in 
Germany.  The  French  army  and  industry  in  France  and  Germany  showed  great  interest  in 
thermoacoustic  refrigeration  technology.  Copies  of  papers  were  distributed  to  scientists  in 
Europe. 

2.  RESEARCH  ACCOMPLISHMENTS 

The  theory,  research  and  applications  of  heat  transfer  equipment,  such  as  thermoacoustic 
refrigerators  or  heat  exchangers  for  thermoacoustic  refrigerators,  involve  a  variety  of  tasks  and 
problems.  In  general,  two  fundamentally  different  tasks  can  be  distinguished  both  from  the 
engineering  as  well  as  theoretical  point  of  view;  namely  the  design  problem  and  the  performance 
prediction.  The  first  task,  the  design  problem,  which  is  frequently  called  the  sizing  problem, 
deals  with  a  design  of  a  new  device  for  a  specific  set  of  input  data  specified  by  operating, 
thermophysical  and  material  parameters  as  well  as  geometry  constraints.  The  second  task,  the 
performance  analysis  or  rating  problem,  involves  the  analysis  and  prediction  of  the  performance 
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of  a  device.  Sometimes,  the  performance  problem  is  treated  as  the  performance  sensitivity 
analysis.  This  term  emphasizes  the  fact  that  the  prediction  also  focuses  on  the  performance  of  a 
system  under  conditions  away  from  the  operating  point  for  which  the  device  was  originally 
designed. 

An  effective  design  of  heat  transfer  equipment  typically  involves  optimization;  the  choice 
of  the  optimization  method  and  the  performance  criteria  depends  on  the  specific  character  of  a 
device,  optimization  function  and  constraints.  Optimization  methods  were  frequently 
implemented  to  design  one  of  the  most  widespread  classes  heat  transfer  equipment,  the  heat 
exchangers.  One  of  the  oldest  optimization  methods  is  the  “Least-Material  Optimization”  was 
described  by  Eckert,  Drake,  1950  {Introduction  to  the  transfer  of  heat  and  mass ,  McGraw-Hill , 
New  York  1950,  pp.  32-34).  Similarly,  “Least-Volume  Optimization”  or  “Optimization  for 
Minimum  Flow  Resistance”  can  be  used,  as  discussed  by  Bejan,  1995  {Convection  heat  transfer, 
John  Wiley  &  Sons,  Inc.,  New  York,  1995,  p.132,  139,  202).  Thermodynamic  optimization 
methods  are  based  on  the  minimization  of  entalphy  exchange  irreversibility  {Bejan,  The  concept 
of  irreversibility  in  heat  exchanger  design:  Counterflow  heat  exchangers  for  gas-gas 
applications ,  J.  Heat  Transfer,  Trans.  ASME,  99,  1977,  pp.  374-380;  Sekulic  and  Herman,  One 
approach  to  irreversibility  minimization  in  compact  crossflow  heat  exchanger  design,  Int.  Com. 
Heat  Mass  Transfer,  13,  1986,  pp.  23-32)  and  on  the  entropy  generation  minimization  {Bejan, 
Entropy  generation  through  heat  and  fluid  flow,  John  Wiley  &  Sons,  New  York,  1982;  Bejan  et 
al.  Thermal  design  and  optimization,  Wiley,  New  York,  1996;  and  Bejan  et  al.,  Thermodynamic 
optimization  of  geometry:  T-  and  Y-  shaped  constructs  of  fluid  streams.  Int.  J.  Therm.  Sci .,  39, 
2000,  pp.  949-960).  There  are  numerous  additional  influences  that  can  affect  the  design,  such  as 
economic  criteria  (initial  or  operating  or  total  costs)  as  well  as  the  limitations  imposed  by  the 
available  technologies  (“Design  for  Manufacturability”  -  see  Iyengar,  Bar-Cohen,  Design  for 
manufacturability  of  SISE  parallel  plate  forced  convection  heat  sinks,  7th  Inter  Society 
Conference  on  Thermal  Phenomena,  Therm  2000,  Las  Vegas,  May  2000,  pp.  141-148). 

However,  the  field  of  thermoacoustic  refrigeration  as  well  as  the  design  of 
ThermoAcoustic  Refrigerators  (TAR)  is  rather  young,  in  particular  when  compared  with  the  heat 
exchanger  tradition,  thus  spectrum  of  optimization  methods  for  TAR  is  quite  limited.  Often, 
TAR  design  is  supported  by  experience  and  intuition,  which  are  the  two  main  sources  for  the 
designer’s  ad  hoc  choices.  A  performance  analysis  in  a  step  following  the  design  process  verified 
an  initiation  estimate  (e.g.,  Hofler,  1986).  Both  steps,  the  design  problem  and  performance 
analysis,  are  typically  repeated  in  an  iterative  cycle  when  optimization  of  system  performance  is 
necessary.  Sometimes  a  new  design  can  be  supported  by  the  knowledge  of  previous  devices  or 
designs,  and  represents  an  incremental  modification  of  a  well-tested  solution  (e.g.  Garrett’s 
approach,  1991,  who  used  an  experimental  experience  by  Hofler,  1986).  Another  optimization 
attempt  where  a  cost  function  of  performance  was  minimized  using  the  simplex  method,  was 
proposed  by  Minner  at  al.  (1995).  Recently,  a  design  optimization  algorithm  was  developed  by 
Wetzel  and  Herman  (1997),  where  the  coefficient  of  performance  was  used  as  a  criterion  of  the 
performance  of  the  thermoacoustic  core. 

The  heat  transfer  research  carried  out  in  the  Heat  Transfer  Lab  of  the  Johns  Hopkins  University 
focused  on  the  following  issues: 
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2.1  DESIGN  OPTIMIZATION  OF  THERMOACOUSTIC  REFRIGERATORS  FOR 
MAXIMUM  COEFFICIENT  OF  PERFORMANCE  (COP)  (Ref.  J1) 

The  coefficient  of  performance  and  coefficient  of  performance  compared  to  Carnot’s  {COP  and 
COPR,  respectively)  play  an  important  role  in  a  design  of  the  TAR.  The  first  major  achievement 
of  our  research  is  the  development  of  the  systematic  design  and  optimization  algorithm  for 
thermoacoustic  refrigerators  (Wetzel  and  Herman,  1997)  that  relies  on  the  first  law  of 
thermodynamics  and  maximizes  the  COP  or  COPR.  This  algorithm  is  based  on  new  scaling 
arguments  that  allowed  significant  simplification  of  the  existing  theory  by  reducing  the  number 
of  design  parameters  from  19  dimensional  parameters  to  10  nondimensional  parameters,  using 
the  simplified  linear  model  of  thermoacoustics,  the  short  stack  boundary  layer  approximation. 
The  methodology  is  based  on  the  first  law  of  thermodynamics.  The  analysis  suggests  a  separate 
optimization  of  the  four  main  modules  of  the  thermoacoustic  refrigerator:  (i)  thermoacoustic  core 
(consisting  of  the  stack  plates  and  the  working  fluid  between  the  stack  plates),  (ii)  resonance 
tube,  (iii)  heat  exchangers  and  the  (iv)  acoustic  driver.  An  upper  limit  on  the  thermoacoustic 
refrigerator’s  performance  is  set  by  the  heat  pumping  capacity  of  the  thermoacoustic  core. 
Calculations  of  the  thermoacoustic  core’s  performance  indicate  that  thermoacoustic  refrigeration 
can  achieve  COPRs  competitive  to  commercially  available  refrigeration  systems,  and  predict  up 
to  40-50%  of  Carnot’s  efficiency  for  the  thermoacoustic  core.  The  design  optimization  algorithm 
introduced  in  Ref.  J1  provides  fast  and  simple  engineering  estimates  that  are  essential  for  the 
development  of  prototypes.  This  is  particularly  important  nowadays,  when  commercial 
applications  are  sought.  A  copy  of  the  paper.  Ref.  J 1 ,  is  attached. 

2.2  DESIGN  OPTIMIZATION  OF  THERMOACOUSTIC  REFRIGERATORS  FOR 
MAXIMUM  COOLING  LOAD  (Ref.  J8) 

Although  use  of  the  COP  and  COPR  as  the  performance  criterion  is  thermodynamically  correct 
and  the  importance  of  these  criteria  obvious,  both  COP  and  COPR  are  relative  parameters,  and 
therefore  have  to  be  used  with  caution.  Recent  studies  focusing  on  the  design  of  a  thermoacoustic 
refrigerator  for  the  thermal  control  of  electronics  have  shown  that  the  design  optimization  can 
require  an  absolute  criterion  (Travnicek  and  Herman,  2002,  Ref.  J8).  A  reason  for  this  is  that  a 
TAR  with  relatively  high  COP  as  well  as  COPR  (both  are  desirable)  can  have  rather  small 
cooling  load,  which  is  unacceptable  in  some  applications.  Hence  the  cooling  load  of  the  TAR 
( Qc )  was  chosen  as  the  criterion  for  optimizing  a  TAR  by  maximizing  the  Qc.  Theoretically,  this 
approach  is  based  on  the  knowledge  that  the  maxima  in  COP  and  Qc  do  not  coincide.  However, 
both  COP  and  Qc  can  be  large  for  a  relatively  wide  range  of  system  parameters,  as  pointed  out 
by  Swift  (1988,  p.  1 171)  for  the  thermoacoustic  engine  (prime  mover)  and  for  the  corresponding 
pair  of  relative  and  absolute  parameters,  i.e.  "engine  efficiency"  and  "extracted  work  flux", 
respectively. 

The  TAR  optimization  methodology  based  on  the  maximization  of  Qc  developed  in  the 
Heat  Transfer  Lab  of  the  Johns  Hopkins  University  by  Z.  Travnicek  and  C.  Herman  appears  to 
be  a  logical  and  adequate  approach  in  a  wide  range  of  applications.  Moreover,  the  requirement  of 
high  Qc  -  values  is  intuitively  and  automatically  taken  into  account  by  the  designer.  A  paper 
describing  this  design  methodology  is  currently  in  preparation  (Ref.  J8). 
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2.3  HOLOGRAPHIC  INTERFEROMETRY  IN  THE  VISUALIZATION  AND 
MEASUREMENT  OF  HIGH-SPPED,  UNSTEADY  TEMPERATURE  FIELDS  IN  THE 
VICINITY  OF  THE  THERMOACOUSTIC  STACK  (Refs.  J2,  J3,  J4) 

Temperature  measurement  by  holographic  interferometry  in  the  presence  of  pressure 
variations  (Refs.  J2,  J3  and  J4) 

In  order  to  quantitatively  reconstruct  the  temperature  distributions  in  the  region  of  the 
thermoacoustic  stack,  it  was  necessary  to  develop  a  new  algorithm  for  the  evaluation  of  the 
interference  fringe  patterns  generated  by  HI  to  obtain  accurate  data  in  the  presence  of  pressure 
variations.  The  algorithms  that  were  previously  described  in  the  literature  do  not  account  for  the 
pressure  variations.  Thus,  by  developing  two  new  evaluation  algorithms,  we  have  expanded  the 
applicability  of  HI  to  temperature  measurements  in  situations  where  pressure  variations 
accompany  the  temperature  oscillations.  Copies  of  the  papers  describing  the  measurement  and 
reconstruction  techniques,  Refs.  J2,  J3  and  J4,  are  attached. 

Improving  the  accuracy  of  temperature  measurements  by  holographic  interferometry  by 
using  temperature  as  tracer  (Ref.  J3) 

We  have  also  shown  that  it  is  possible  to  enhance  the  spatial  resolution  and  accuracy  of 
interferometric  measurements  by  imposing  a  steady-state  temperature  distribution  to  the 
oscillating  flow.  The  latter  technique  can  be  considered  as  using  temperature  as  tracer  to  resolve 
the  small,  time  dependent  temperature  fluctuations  in  the  acoustic  field.  A  copy  of  the  paper. 
Ref.  J3,  is  attached. 

2.4  EXPERIMENTAL  STUDY  OF  THERMOACOUSTIC  EFFECTS  ON  A  SINGLE 
STACK  PLATE:  TEMPERATURE  FIELDS  AND  HEAT  TRANSFER  (Refs.  J5  and  J6) 

Experiments  with  the  modular  thermoacoustic  refrigerator 

Experiments  were  carried  out  with  a  modular  thermoacoustic  refrigerator  to  study  the  thermal 
performance  of  the  stack  region.  To  investigate  the  influence  of  key  parameters  of  the 
thermoacoustic  refrigerator  on  heat  transfer,  we  built  a  modular  thermoacoustic  refrigerator 
based  on  the  design  algorithm  developed  by  Wetzel  and  Herman  (1997-J1)  which  represents  the 
design  parameters  in  dimensionless  form.  Introducing  dimensionless  design  parameters  has  the 
advantage  that  the  results  are  not  restricted  to  one  particular  experimental/hardware 
configuration.  Additionally,  the  approach  allows  comparison  with  existing  designs  and  reduces 
the  number  of  experiments  that  need  to  be  conducted.  The  modular  structure  of  the  refrigerator 
allows  easy  modification  of  stack  and  heat  exchanger  geometries,  such  that  the  influence  of  the 
different  system  parameters  on  the  heat  transfer  can  be  investigated. 

In  addition  to  conventional  sensors,  HI  combined  with  high-speed  cinematography  was 
applied  as  measurement  technique,  as  described  in  Section  2.3.  The  main  advantage  of  this 
approach  is  that  it  allows  high  spatial  (up  to  2700  dpi)  and  high  temporal  resolutions  (up  to 
0.1ms).  Thus  it  represents  the  ideal  tool  to  investigate  the  temperature  fields  in  a  thermoacoustic 
refrigerator  that  oscillate  with  frequencies  of  a  few  hundred  hertz. 

We  designed  and  built  three  stacks  for  the  modular  thermoacoustic  refrigerator.  Stack  1 
was  used  for  temperature  measurements  with  thermocouples  as  well  as  HI.  Stack  2  was  identical 
in  its  dimensions  and  instrumentation  to  stack  1,  but  had  copper  strips  that  act  as  heat  exchangers 
attached  at  both  ends.  Stack  3  was  also  identical  in  its  dimensions  to  stack  1,  but  instead  of 


thermocouples  it  was  equipped  with  a  pressure  transducer.  The  pressure  transducer  measures 
pressure  at  the  stack  center  position.  All  three  stacks  were  movable  within  the  resonance  tube, 
and  therefore  temperature  and  pressure  measurements  for  different  stack  center  positions  could 
be  performed.  By  conducting  temperature  measurements  with  stacks  1  and  2,  the  effect  of  the 
heat  exchangers  on  the  temperature  distributions  was  investigated.  Pressure  measurements  with 
stack  3  allow  quantifying  the  effect  of  the  stack  on  the  acoustic  standing  wave.  The  pressure 
distribution  along  the  acoustic  standing  wave  can  be  determined  by  measuring  the  pressure  at 
different  locations  along  the  resonance  tube  with  a  second  pressure  transducer. 

Visualization  and  measurement  of  the  oscillating  temperature  fields  at  the  edge  of  the 
thermoacoustic  stack  plates  (Refs.  J5  and  J6) 

The  second  major  achievement  of  the  experimental  part  of  the  project  is  the  measurement  of 
oscillating  temperature  fields  as  well  as  heat  fluxes  in  a  thermoacoustic  refrigerator  model  using 
a  combination  of  real-time  holographic  interferometry  (HI)  and  high-speed  cinematography.  The 
thermal  interaction  between  a  heated  solid  plate  and  the  acoustically  driven  working  fluid  was 
investigated  in  the  Heat  Transfer  Laboratory  of  the  Johns  Hopkins  University  by  visualizing  and 
quantifying  the  temperature  fields  in  the  neighborhood  of  the  solid  plate.  With  these  experiments 
the  thermoacoustic  effect  was  visualized  for  the  first  time.  The  measurements  the  existence  of 
high  time-averaged  local  heat  fluxes  at  the  edge  of  the  stack  plates,  as  predicted  by  numerical 
and  theoretical  treatments  of  the  problem.  A  better  knowledge  of  these  temperature  fields  is 
essential  to  develop  systematic  design  methodologies  for  heat  exchangers  in  oscillatory  flows. 

Key  results  of  these  studies  were  reported  in  two  papers,  refs.  J5  and  J6  (attached).  The 
difference  between  heat  transfer  in  oscillatory  flows  with  zero  mean  velocity  and  steady  flows  is 
discussed  in  the  papers.  In  the  experiments,  the  thermoacoustic  effect  was  visualized  through 
temperature  measurements.  A  novel  evaluation  procedure  that  accounts  for  the  influence  of  the 
acoustic  pressure  variations  on  the  refractive  index  was  applied  to  accurately  reconstruct  the 
high-speed,  two-dimensional  oscillating  temperature  distributions.  Heat  fluxes  close  to  the  edge 
of  a  heated  solid  plate  aligned  parallel  to  the  axis  of  the  acoustic  standing  wave  were  measured 
for  drive  ratios  DR  =  PA/pm  of  1,  2  and  3%.  At  the  highest  drive  ratio  (3%),  the  resulting  heat 
flux  vector  at  the  edge  of  the  plate  is  directed  into  the  plate,  opposite  to  the  direction  of  the  heat 
flux  imposed  by  the  resistive  heaters  within  the  plate.  This  observation  confirms  the 
thermoacoustic  effect  previously  discovered  in  the  visualized  temperature  fields.  Through  the 
energy  balance  the  magnitudes  of  the  heat  fluxes  into  the  plate,  caused  by  the  thermoacoustic 
effect,  were  determined.  The  measured  data  are  in  good  agreement  with  numerical  and 
analytical  predictions. 

2.5  IMPACT  OF  FLOW  OSCILLATIONS  ON  CONVECTIVE  HEAT  TRANSFER  (Ref. 

B2) 

Flow  oscillations  and  their  impact  on  heat  transfer  from  solid  surfaces  play  a  key  role  in  the 
modeling  and  performance  of  heat  exchangers  for  thermoacoustic  applications.  A  review  paper 
in  the  form  of  a  book  chapter,  discussing  the  impact  of  flow  oscillations  on  heat  transfer  in  self- 
sustained  oscillatory  flows,  jets  and  some  acoustically  driven  flows  was  prepared  (Ref.  B2, 
attached). 

2.6  A  SIMPLIFIED  MODEL  OF  HEAT  TRANSFER  IN  PARALLEL  CHANNEL  HEAT 
EXCHANGERS  AND  STACK  PLATES  OF  THERMOACOUSTIC  DEVICES  (Ref.  J7) 


The  heat  exchangers  of  a  thermoacoustic  refrigerator  deliver  the  heat  absorbed  in  the  refrigerated 
volume  to  the  cold  side  of  the  thermoacoustic  stack  and  absorb  the  heat  transported  from  the  cold 
to  the  hot  side  of  the  stack,  which  is  then  released  into  the  ambient  To  improve  the  performance 
of  the  heat  exchangers,  it  is  essential  to  gain  better  insight  into  the  thermoacoustic  heat  transfer 
process.  So  far,  neither  theoretical  nor  numerical  models  are  available  that  can  accurately  predict 
the  overall  heat  transfer  performance  of  heat  exchangers  as  part  of  thermoacoustic  devices 
(numerical  results  regarding  temperature  and  velocity  distributions  around  the  stack  plates  were 
reported  in  the  literature)  and  their  influence  on  the  overall  performance  of  the  TAR.  Therefore, 
one  of  the  goals  of  the  project  involved  quantifying  thermoacoustic  heat  transfer  in  the  heat 
exchangers  in  a  combined  modeling  and  experimental  effort. 

Standard  heat  exchanger  design  methods  are  available  for  conventional  applications 
(steady  forced  convection)  for  a  wide  range  of  heat  exchanger  designs,  flow  arrangements  and 
working  fluids.  Up  to  date  it  has  not  been  evaluated  to  which  extent  conventional  methods  are 
suitable  for  the  design  of  heat  exchangers  operating  in  oscillatory  flows  with  zero  mean  velocity. 
Currently  it  is  being  explored  if  it  is  feasible  to  modify  existing  methods  as  opposed  to 
developing  new  models  and  methods  suitable  for  predicting  the  behavior  thermoacoustic  heat 
exchangers.  The  theoretical  analysis  will  provide  information  on  the  experimental  data, 
regarding  both  heat  transfer  and  pressure  drop,  that  are  necessary  for  the  design  and  modeling  of 
the  heat  exchangers.  Finally,  the  heat  exchanger  models  will  be  incorporated  into  the 
optimization  algorithm  described  in  Jl. 

Mathematical  models  of  the  thermal  behavior  of  several  types  of  heat  exchangers  were 
developed  for:  (i)  the  parallel  flat-channel  heat  exchanger  using  a  liquid  as  the  transport  fluid 
(J7),  (ii)  the  simple,  parallel-fin  heat  exchanger  that  relies  on  heat  conduction  as  well  as  (iii)  2D 
computational  models  for  several,  more  complex  heat  exchanger  configurations. 

A  simplified  model  of  heat  transfer  in  parallel  channel  heat  exchangers  and  stack  plates 
of  thermoacoustic  devices  was  developed.  The  model  took  advantage  of  previous  results 
regarding  the  thermal  behavior  of  the  thermoacoustic  core  and  the  measurements  of  local  heat 
transfer  at  the  edge  of  the  stack  plates  for  investigations  of  the  performance  of  heat  exchangers 
attached  to  the  thermoacoustic  core.  Geometrical  and  operational  parameters  as  well  as 
thermophysical  properties  of  the  heat  exchangers,  the  stack  plate,  and  the  working  fluid  were 
organized  into  dimensionless  groups  that  allowed  accounting  for  their  impact  on  the  performance 
of  the  heat  exchangers.  Numerical  simulations  were  carried  out  with  the  model.  Nonlinear 
temperature  distributions  and  heat  fluxes  were  observed  near  the  edge  of  the  stack  plate.  Effects 
of  different  parameters  on  the  thermal  performance  of  the  heat  exchangers  were  investigated. 

2.7  DEVELOPMENT  OF  A  SIMPLIFIED  MODEL  OF  THE  HEAT  EXCHANGERS 
RELYING  ON  HEAT  CONDUCTION  -  PARALLEL  STRIP  HEAT  EXCHANGERS  (Ref. 
J9) 

The  simple,  parallel  fin  heat  exchanger  design  that  relies  on  heat  conduction  to  transport  heat  to 
and  from  the  thermoacoustic  core  delivers  satisfactory  heat  transfer  rates  for  smaller 
thermoacoustic  refrigerators  because  of  the  small  cross  sectional  area  of  the  resonance  tube.  The 
mechanism  of  heat  transfer  along  the  fins  is  conduction.  The  conductive  resistances  of  the  fins 
scale  with  the  fin  length,  and  are  small  for  miniature  systems.  This  solution  for  the  design  of  the 
heat  exchangers  offers  the  advantage  of  simplicity,  it  eliminates  the  problem  of  sealing  at  the 
locations  where  the  pipe  with  the  heat  exchanger  liquid  enters  the  resonance  tube  and  the 
manufacturing  of  the  heat  exchanger  hardware  is  easier  than  for  heat  exchangers  that  use  liquids. 


In  the  first  step  we  assume  that  the  temperature  of  the  heat  exchanger  fin  is  uniform  along 
the  acoustic  axis  x.  This  assumption  is  justifiable:  the  thermal  conductivity  of  the  fin  material 
(typically  copper)  is  high  and  the  dimension  along  the  acoustic  axis  is  small.  The  steady  state 
temperature  distribution  along  the  fin  was  determined  analytically  and  numerically.  Heat  transfer 
from  the  oscillating  flow  (with  zero  mean  velocity)  is  described  by  means  of  an  equivalent 
constant  heat  transfer  coefficient  h  (measured,  numerically  determined  or  estimated  values  were 
used  in  the  calculations).  The  ring  that  spans  the  heat  exchanger  fins  is  assumed  to  be  at  the 
ambient  temperature,  which  corresponds  to  the  temperature  of  the  fin  base.  The  solution  of  the 
model  yields  temperature  and  heat  flux  distributions  in  the  heat  exchangers.  Initially  a  solution 
for  a  single  fin  (the  longest  one)  was  determined;  temperature  distributions  in  the  shorter  fins 
were  found  using  appropriate  scaling. 

2.8  DESIGN  OF  A  MINIATURE  THERMOACOUSTIC  REFRIGERATOR  FOR 
ELECTRONIC  COOLING  APPLICATIONS  (Ref.  J9) 

Design  issues  were  considered  that  influence  the  performance  of  miniature  thermoacoustic 
refrigerators  with  application  to  the  cooling  of  electronic  equipment.  In  the  design  optimization 
process  the  design  objective  was  to  maximize  the  cooling  load.  It  was  found  that  the  design  that 
maximizes  the  cooling  load  is  not  identical  to  the  design  that  maximizes  the  coefficient  of 
performance  of  the  system,  COP.  It  was  found  that  Helium  as  working  fluid  will  have  the  highest 
cooling  load.  The  influence  of  the  geometry  and  thermophysical  parameters  on  the  performance 
of  the  thermoacoustic  refrigerator  was  studied  systematically  and  a  performance  sensitivity 
analysis  was  carried  out.  Key  results  are  being  summarized  in  a  paper  by  Travnicek  and  Herman 
(2002),  Ref.  J9. 

2.9  INCORPORATION  OF  HEAT  EXCHANGER  MODELS  INTO  THE  DESIGN 
OPTIMIZATION  ALGORITHM  BY  WETZEL  AND  HERMAN  (1997) 

The  models  of  the  heat  exchangers  will  be  incorporated  into  the  design  optimization  algorithms 
developed  by  Wetzel  and  Heiman  (1997),  Ref.  Jl,  and  Travnicek  and  Herman  (2001),  Ref.  J9,  to 
allow  the  quantification  of  the  impact  of  heat  exchangers  on  the  overall  performance  of 
thermoacoustic  devices.  Performance  limitations  of  miniature  thermoacoustic  refrigerators  will 
be  evaluated. 

2.10  FUTURE  WORK 

Measurement  of  temperature  and  velocity  profiles  within  the  viscous  and  thermal 
penetration  depths 

One  major  research  issue  that  needs  to  be  addressed  is  to  experimentally  determine  the  time 
dependent  temperature  and  velocity  profiles  within  the  thermal  and  viscous  penetration  depths, 
respectively.  This  knowledge  is  important  for  two  reasons.  First,  practical  thermoacoustic 
devices  will  most  likely  operate  in  non-linear  regimes.  Therefore,  it  is  necessary  to  identify  the 
influence  of  non-linear  effects  on  the  phasing  between  temperature  and  velocity  within  the 
thermal  penetration  depth,  since  it  is  this  phasing  that  generates  the  hydrodynamic  transport  of 
heat  along  the  stack  plates.  Second,  in  order  to  accurately  measure  the  time  dependent  heat 
fluxes  within  the  stack  and  the  heat  exchangers  in  thermoacoustic  devices,  it  is  essential  to 
accurately  determine  the  derivative  of  the  temperature  profile  near  the  solid  boundary.  It  is 
planned  to  acquire  a  Laser  Doppler  Anemometer  to  measure  the  velocity  profile  within  the 
viscous  penetration  depth. 


The  data  obtained  in  these  experiments  will  be  compared  quantitatively  to  results  of 
numerical  simulations  (obtained  by  Omar  Knio,  for  example)  to  verify  the  numerical  model.  It  is 
expected  that  the  numerical  model  will  provide  fundamental  data  necessary  for  the  optimization 
of  heat  exchangers.  However,  comparisons  with  experimental  data  will  be  necessary  after  any 
significant  change  in  the  numerical  algorithm  or  physical  model. 

Development  of  a  micro  heat  pipe  serving  as  heat  exchanger  in  thermoacoustic 
refrigerators 

Apart  from  conventional  heat  exchangers  heat  pipes  can  also  be  implemented  to  remove  heat 
from  the  refrigerated  volume  and  deliver  it  to  the  thermoacoustic  refrigerator.  Heat  pipes  offer 
the  advantage  of  being  passive  devices,  which  makes  them  particularly  suitable  for  applications 
where  high  reliability  is  essential.  A  micro-heat  pipe  serving  as  heat  exchanger  will  be  designed 
and  built  based  on  the  information  regarding  pressure,  velocity  and  temperature  distributions 
obtained  experimentally. 

Improving  the  design  optimization  algorithms  and  incorporating  heat  exchanger  models 
into  the  performance  calculation 

The  optimization  criteria  used  in  the  development  and  design  of  thermoacoustic  devices  depend 
on  the  application,  as  described  in  sections  2.1-2.2  (Refs.  J1  and  J8).  Additional  criteria  may  be 
required  in  novel  applications  as  well  as  when  coupling  thermoacoustic  prime  movers  with 
thermoacoustic,  Stirling  or  pulse  tube  refrigeration  systems.  A  systematic  investigation  and 
optimization  of  the  coupling  requires  further  study.  Also,  incorporating  models  of  the  heat 
exchangers  into  the  design  optimization  algorithms  will  lead  to  further  enhancement  of  the 
performance  of  thermoacoustic  systems,  which  is  essential  in  practical  implementations. 
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3.1  Peer-reviewed  journal  papers 

Jl.  Wetzel,  M.,  Herman,  C.,  1997,  Design  optimization  of  thermoacoustic  refrigerators , 
International  Journal  of  Refrigeration,  Vol.  20,  No.l,  pp.  3-21 . 

J2.  Wetzel,  M.,  Herman,  C.,  1998,  Accurate  measurement  of  high-speed,  unsteady  temperature 
fields  by  holographic  interferometiy  in  the  presence  of  periodic  pressure  variations , 
Measurement  Science  and  Technology,  Vol.  9,  No.  6,  pp.  939-951. 

J3.  Herman,  C.,  Kang,  E.,  Wetzel,  M.,  1998,  Expanding  the  applications  of  holographic 
interferometry  to  the  quantitative  visualization  of  complex ,  oscillatory  thermofluid  processes , 
Experiments  in  Fluids,  Vol.  24,  pp.  431-446. 

J4.  Wetzel,  M.,  Herman,  C.,  1998,  Limitations  of  temperature  measurements  with  holographic 
interferometry  in  the  presence  of  pressure  variations ,  Experimental  Thermal  and  Fluid 
Science,  Vol.  17,  pp.  294-308. 

J5.  Wetzel,  M.,  Herman,  C.,  2000,  Experimental  study  of  thermoacoustic  effects  on  a  single 
stack  plate  -  Parti:  Temperature  fields ,  Heat  and  Mass  Transfer,  Vol.  36,  pp.  7-20. 

J6.  Wetzel,  M.,  Herman,  C.,  1999,  Experimental  study  of  thermoacoustic  effects  on  a  single 
stack  plate  -  Part  II:  Heat  Transfer ,  Heat  and  Mass  Transfer,  Vol.  35,  pp.  433-441 . 

J7.  Chen,  Y.,  Herman,  C.,  2001,  A  simplified  model  of  heat  transfer  in  parallel,  flat-channel  heat 
exchangers  for  thermoacoustic  applications,  accepted  for  publication  in  Heat  and  Mass 
Transfer. 


J8.  Travnicek,  Z.,  Herman,  C.,  2002,  Optimization  of  thermoacoustic  refrigerators  for  maximum 
cooling  load ,  in  preparation. 

19.  Travnicek,  Z.,  Herman,  C.,  2002,  Thermoacoustic  refrigerator  for  the  cooling  of  electronic 
equipment ,  in  preparation. 

3.2  Invited  and  peer-reviewed  book  chapters 

Bl.  Herman,  C.,  1999,  Heat  exchangers  for  thermoacoustic  refrigerators:  Heat  transfer  issues  in 
oscillatoiy  flow,  Energy  Conservation  through  Heat  Transfer  Enhancement  of  Heat 
Exchangers,  Editors:  S.  Kakac  et  al.,  NATO  Advanced  Science  Institute  Series  E,  Vol.  355, 
Kluwer  Academic  Publishers,  Dodrecht,  pp.  285-299,  peer-reviewed  and  invited 
contribution. 

B2.  Herman,  C.,  2000,  The  impact  of  flow  oscillations  on  convective  heat  transfer ;  Annual 
Review  of  Heat  Transfer,  Editor:  C.-L.  Tien,  Vol.  XI,  Chapter  8,  pp.  495-562,  invited 
contribution. 

B3.  Herman,  C.,  2001,  one  section  in  the  Dictionary  of  Acoustics,  Editor:  D.  Basu,  CRC  Press, 
Boca  Raton,  invited  contribution. 

3.3  Conference  papers 

Cl. Herman,  C.,  Wetzel,  M.,  Leungki  Y.  S.,  1995,  Visualization  of  oscillating  temperature  and 
flow  fields  in  a  thermoacoustic  refrigerator ,  Proceedings  of  the  Seventh  International 
Symposium  on  Flow  Visualization,  Seattle  Washington,  Begell  House,  Inc.  New  York,  pp. 
334-340. 

C2.  Herman,  C.,  Wetzel,  M.,  1995,  Design  of  a  thermoacoustic  refi'igerator  -  a  case  study-,  AES- 
Vol.35,  ASME,  pp.  195-203. 

C3.  Wetzel,  M.,  Herman,  C.,  1996,  Design  issues  of  a  thermoacoustic  refrigerator  and  its  heat 
exchangers,  ASME  Proc.  of  the  31st  National  Heat  Transfer  Conference,  HTD-Vol.  331,  pp. 
137-144. 

C4.  Wetzel,  M.,  Herman,  C.,  1996,  Parameter  spaces  and  design  optimization  of  thermoacoustic 
refrigerators ,  AES-Vol.  36:  1071-6947,  Proc.  1996  ASME  IMECE,  Symposium  on 
Thermodynamics  and  the  Design,  Analysis  and  Improvement  of  Energy  Systems,  Atalanta, 
pp.  355-363. 

ASME  Student  Design  Paper  Award. 

C5.  Wetzel,  M.,  C.  Herman,  Accurate  measurements  of  temperature  fields  by  holographic 
interferometry  in  the  presence  of  pressure  variations,  Proc.  8th  Int.  Symp.  on  Flow 
Visualization,  Editor  Ian  Grant,  Sorrento,  Italy,  Sept.  1998. 

C6.Chen,  Y.,  C.  Herman,  1999,  A  simplified  model  of  heat  transfer  in  heat  exchangers  and  stack 
plates  of  thermoacoustic  devices,  Proc.  National  Heat  Transfer  Conference,  Albuquerque, 
NM,  Paper  NHTC99-3 10,  pp.  1-8. 

3.4  Conference  proceedings  -  abstract  published 

C7. Herman,  C.,  Bartscher,  C.,  Wetzel,  M.,  Volejnik,  M.,  1994,  Experimental  visualization  of 
heat  transfer  and  fluid  flow  processes  in  a  thermoacoustic  device,  J.  Acoust.  Soc.  Am,  Vol. 
96,  No.  5,  Pt.  2,  pp.  3220-3221. 

C8. Herman,  C.,  Wetzel,  M.,  Volejnik,  M.,  1995,  Visualization  of  oscillating  temperature  and 
flow  fields  in  the  stack  region  of  a  thermoacoustic  refrigerator  model,  J.  Acoust.  Soc.  Am, 
Vol.  97,  No.  5,  Pt.  2,  p.  3410. 


C9.Herman,  C.,  Wetzel,  M.,  Wagner,  J.,1995,  Holographic  interferometry:  An  approach  to  study 
the  unsteady  temperature  fields  in  the  stack  region  and  its  neighborhood,  J.  Acoust.  Soc. 
Am,  Vol.  98,  No.  5,  Pt.  2,  p.  2962. 

CIO.  Wetzel,  M.,  Herman,  C.,1996,  Optimization  of  the  performance  of  thermoacoustic 
refrigerators  applying  the  short  stack  boundary  layer  approximation ,  J.  Acoust.  Soc.  Am, 
Vol.  99,  No.  4,  Pt.  2,  p.  2559. 

Cll.  Wetzel,  M.,  Herman,  C.,  1996,  Measurements  of  oscillating  temperature  fields  in  the 
stack  region  of  a  thermoacoustic  refrigerator  model ,  J.  Acoust.  Soc.  Am,  Vol.  100,  No.  4,  Pt. 
2,  p.  2846. 

C12.  Wetzel,  M.,  Herman,  C.,  1997,  Heat  transfer  measurements  in  a  thermoacoustic 
refrigerator ,  J.  Acoust  Soc.  Am,  Vol.  101,  No.  5,  Pt.  2,  pp.  3021-3022. 

C13.  Wetzel,  M.,  Herman,  C.,  1997,  Design  of  a  thermoacoustic  refrigerator  for  visualization 
measurements,  J.  Acoust.  Soc.  Am,  Vol.  102,  No.  5,  Pt.  2,  pp.  3071-3072. 

C14.  Wetzel,  M.,  Herman,  C.,  1998,  Heat  transfer  measurements  on  a  single-plate 
thermoacoustic  refrigerator ;  136th  Meeting  of  the  Acoustical  Society  of  America,  Norfolk, 
VA,  Oct.  1998,  J.  Acoust.  Soc.  Am.,  Vol.  104,  No.  3,  Pt.  2,  Sept.  1998. 

Cl 5.  Herman,  C.,  Conventional  heat  exchanger  design  methods  and  their  applicability  to 
thermoacoustics ,  Acta  Acoutica,  Vol.  85,  Suppl.  1,  p.  86,  1999. 

3.5  Dissertation 

Dl.M.  Wetzel,  Experimental  investigation  of  heat  transfer  in  a  single-plate  thermoacoustic 
refrigerator,  thesis  defense  scheduled  July  14, 1998. 

4.  LECTURES  AND  PRESENTATIONS  AT  CONFERENCES 

1.  Experimental  visualization  of  heat  transfer  and  fluid  flow  processes  in  a  thermoacoustic 
device ,  Herman,  C.,  Bartscher,  C.,  Wetzel,  M.,  Volejnik,  128th  Meeting  of  Acoust.  Soc.  Am, 
Austin,  TX,  Nov.  28-Dec.  2,  1994. 

2.  Visualization  of  oscillating  temperature  and  flow  fields  in  the  stack  region  of  a 
thermoacoustic  refrigerator  model,  Herman,  C.,  Wetzel,  M.,  Volejnik,  M.,  129th  Meeting  of 
Acoust.  Soc.  Am,  Washington,  D.C.,  May  30- June  03,  1995. 

3.  Visualization  of  oscillating  temperature  and  flow  fields  in  a  thermoacoustic  refrigerator , 
Herman,  C.,  Wetzel,  M.,  Leungki  Y.  S.,  Seventh  International  Symposium  on  Flow 
Visualization,  Seattle  Washington,  1995. 

4.  Design  of  a  thermoacoustic  refrigerator  -  a  case  study,  Herman,  C.,  Wetzel,  M.,  IMECE,  San 
Francisco,  November  1995. 

5.  Holographic  inteiferometiy:  An  approach  to  study  the  unsteady  temperature  fields  in  the 
stack  region  and  its  neighborhood,  Herman,  C.,  Wetzel,  M.,  Wagner,  J,  130th  Meeting  of 
Acoust.  Soc.  Am,  St.  Louis,  MO,  Nov.  27-Dec.01,  1995. 

6.  Visualization  of  high-speed  heat  transfer  phenomena,  Department  of  Mechanical 
Engineering,  University  of  Maryland  Baltimore  County,  Baltimore,  December  8,  1995. 

7.  A  method  of  quantitative  visualization  of  complex,  unsteady  thermofluid  processes,  Herman, 
C.,  Kang,  E.,  ASME  Fluids  Engineering  Division  Conference,  San  Diego,  1996. 

8.  Design  issues  of  a  thermoacoustic  refrigerator  and  its  heat  exchangers,  Wetzel,  M.,  Herman, 
C.,  3 1st  National  Heat  Transfer  Conference,  Houston,  TX,  August  03-06,  1996L 

9.  Parameter  spaces  and  design  optimization  of  thermo  acoustic  refrigerators,  Wetzel,  M., 
Herman,  C.,  IMECE,  Atalanta,  November  1996. 


10.  Optimization  of  the  performance  of  thermoacoustic  refrigerators  applying  the  short  stack 
boundary  layer  approximation ,  Wetzel,  M.,  Herman,  C.,  131th  Meeting  of  Acoust.  Soc.  Am, 
Indianapolis,  IN,  May.  13-17,  1996. 

1 1 .  Measurements  of  oscillating  temperature  fields  in  the  stack  region  of  a  thermoacoustic 
refrigerator  model ,  Wetzel,  M.,  Herman,  C.,  132th  Meeting  of  Acoust.  Soc.  Am,  Honolulu, 
HI.  Dec.  02-06,1996. 

12.  Heat  transfer  measurements  in  a  thermoacoustic  refrigerator,  Wetzel,  M.,  Herman,  C.,  133th 
Meeting  of  Acoust.  Soc.  Am,  State  College,  PA,  June  15-20,  1997. 

13.  Design  of  a  thermoacoustic  refrigerator  for  visualization  measurements,  Wetzel,  M., 
Herman,  C.,  134th  Meeting  of  Acoust.  Soc.  Am,  San  Diego,  CA.  Dec.  1,1997. 

14.  Simultaneous  measurement  and  visualization  of  oscillatory  temperature  and  flow  fields  by 
holographic  interferometry ,  Ecole  Politechnique  Federale  de  Lausanne,  Lausanne, 
Switzerland,  January  12,  1998. 

1 5.  Advances  in  the  study  of  forced  convection  heat  transfer  processes  using  holographic 
interferometry,  National  Institute  of  Standards  and  Technology,  Gaithersburg,  MD,  Feb.  3, 

1998. 

16.  Cool  sound  -  The  history  and  future  of  thermoacoustic  refrigeration,  Baltimore  Chapter  of 
ASHRAE,  Feb.  19, 1998. 

17.  Using  temperature  as  tracer:  new  possibilities  for  an  “old”  tool  -  real-time  holographic 
interferometry  -  in  the  study  of  high-speed  thermofluid  processes,  University  of  Maryland, 
College  Park,  Feb.  27, 1998. 

18.  Thermoacoustic  refrigeration,  Women’s  Science  Forum,  Space  Telescope  Science  Institute, 
Baltimore,  Feb.  28, 1998. 

19.  Experimental  investigation  of  thermoacoustic  refrigeration,  Colorado  State  University,  Fort 
Collins,  March  5,  1998. 

20.  Experimental  investigation  of  thermoacoustic  refrigeration,  University  of  Florida, 
Gainesville,  March  20, 1998. 

21.  Heat  exchangers  for  thermoacoustic  refrigerators :  Heat  transfer  issues  in  oscillatory  flow, 
invited  lecture,  NATO  Advanced  Study  Institute  series  on  Energy  Conservation  through  Heat 
Transfer  Enhancement  of  Heat  Exchangers,  May  25  -  June  5,  1998,  Cesme/Izmir,  Turkey. 

22.  Quantitative  visualization  of  high-speed  convective  heat  transfer  processes,  Dynaflow 
Corporation,  Baltimore,  Sept.  1998. 

23.  Heat  transfer  measurements  on  a  single-plate  thermoacoustic  refrigerator,  136th  Meeting  of 
the  Acoustical  Society  of  America,  Norfolk,  VA,  Oct.  1998. 

24.  Quantitative  visualization  of  high-speed  heat  transfer  phenomena:  Applications  in 
thermoacoustics.  The  Jamie  Whitten  National  Center  for  Physical  Acoustics,  University  of 
Mississippi,  Dec.  8, 1998. 

25.  Heat  exchangers,  conventional  methods  of  heat  exchanger  design ,  and  their  role  and 
applicability  in  thermoacoustics,  The  Jamie  Whitten  National  Center  for  Physical  Acoustics, 
University  of  Mississippi,  Dec.  10,  1998. 

26.  Experimental  investigation  of  thermoacoustic  refrigeration,  University  of  California  at  Los 
Angeles,  Los  Angeles,  December  14,  1998. 

27.  Progress  in  the  experimental  investigation  of  thermoacoustic  refrigeration,  invited  lecture, 
Thermoacoustics  Review  Meeting  —  TARM  —  organized  by  the  Office  of  Naval  Research, 
The  Jamie  Whitten  National  Center  for  Physical  Acoustics,  University  of  Mississippi,  Jan.  7, 

1999. 
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28.  Experimental  investigation  of  thermoacoustic  refrigeration ,  Drexel  University,  Philadelphia, 
January  29,  1999. 

29.  A  simplified  model  of  heat  transfer  in  heat  exchangers  and  stack  plates  of  thermoacoustic 
refrigerator ,  Joint  Meeting:  ASA/EAA/DEGA,  Berlin,  March  16,  1999. 

30.  Thermoacoustic  refrigerators:  Design  optimization  and  experimental  investigation ,  Reunion 
ThermoAcoustique,  Meeting  of  the  French  Thermoacoustics  Research  group,  Ecole  Centrale 
de  Lyon,  France,  March  24,  1999. 

31  .Experimental  investigation  of  theimoacoustic  refrigeration ,  Ecole  Centrale  Paris,  France, 
March  26, 1999. 

32.  Theimoacoustic  refrigeration:  Design  optimization ,  experiments  and  miniaturization , 
Rockwell  Science  Corporation,  CA,  April  12,  1999. 

33.  Experimental  investigation  of  thermoacoustic  refrigeration ,  NASA  Glenn  Research  Center, 
Cleveland,  OH,  Sept.  8,  1999. 

34.  Experimental  investigation  of  thermoacoustic  refrigeration,  Case  Western  Reserve 
University,  Cleveland,  OH,  Sept.  10,  1999. 

35.  Design  optimization  of  thermoacoustic  refrigerators ,  Universite  Pierre  et  Marie  Curie  (Paris 
6)  -  CNRS  -  LIMSI,  Orsay,  France,  June  7,  2000. 

36.  Expanding  the  application  of  holographic  interferometry  to  the  measurement  of  high-speed, 
unsteady  temperature  distributions  In  the  presence  of  pressure  variations  -  application  to 
thermoacoustics,  Universite  Pierre  et  Marie  Curie  (Paris  6)  —  CNRS  -  LIMSI,  Orsay,  France, 
June  14,  2000. 

37.  Experimental  study  of  thermoacoustic  effects  on  a  single  stack  plate:  temperature  fields  and 
heat  transfer ;  Universite  Pierre  et  Marie  Curie  (Paris  6)  —  CNRS  -  LIMSI,  Orsay,  France, 
June  21, 2000. 

38.  A  simplified  model  of  heat  transfer  in  parallel  channel  heat  exchangers  and  stack  plates  of 
thermoacoustic  devices,  Universite  Pierre  et  Marie  Curie  (Paris  6)  —  CNRS  -  LIMSI,  Orsay, 
France,  June  23,  2000. 

39.  Performance  of  miniature  thermoacoustic  refrigerators,  Universite  Pierre  et  Marie  Curie 
(Paris  6)  -  CNRS  -  LIMSI,  Orsay,  France,  June  26,  2000. 

40.  Visualization  of  high-speed  flow  and  temperature  fields  in  the  stack  of  a  thermoacoustic 
refrigerator,  National  Heat  Transfer  Conference  Pittsburgh,  August  20,  2000. 

4 1 .  Advances  in  the  investigation  of  thermoacoustic  refrigeration:  modeling  and  experiments, 
University  of  Michigan,  Ann  Arbor,  January  19,  2001. 

42.  Verbesserung  thermoakustischer  Kaltemaschinen  (Improving  the  performance  of 
thermoacoustic  refrigerators ),  Thermodynamics  Colloquium,  University  of  Hannover, 
Hannover,  Germany,  Nov.  23,  2001. 

43.  Thermoacoustic  refrigeration  for  electronics  cooling,  Navy  Electronics  Cooling  Workshop, 
US  Naval  Academy,  Annapolis,  MD,  Feb.  26,  2002. 

44.  Low-temperature  applications  of  thermoacoustic  refrigeration,  invited  lecture  to  be 
presented  at  the  NATO  Advanced  Study  Institute  on  "Low-temperature  and  cryogenic 
refrigeration  -  fundamentals  and  applications",  June  23  -  July  5,  2002,  Altin  Yunus-^esme, 
Izmir,  Turkey. 

45.  Refrigeration  systems  for  low -temperature  electronics,  invited  lecture  to  be  presented  at  the 
NATO  Advanced  Study  Institute  on  "Low-temperature  and  cryogenic  refrigeration  - 
fundamentals  and  applications",  June  23  -  July  5,  2002,  Altin  Yunus-Qesme,  Izmir,  Turkey. 
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5.  AWARDS  AND  HONORS  DURING  PROJECT  PERIOD 

During  the  project  period  both  the  co-PI,  Cila  Herman,  and  the  graduate  student  employed  on 

this  project,  Martin  Wetzel,  received  awards  which  are  listed  below. 

5.1  Awards:  Cila  Herman 

Major  awards 

1.  Presidential  Early  Career  Award  for  Scientists  and  Engineers:  PECASE  -  October  1997, 
Nominating  agency:  NASA.  For  details  see  under  Grants  and  Contracts. 

2.  National  Science  Foundation  (NSF)  -  CAREER  Award  -  July  1997.  For  details  see  under 
Grants  and  Contracts. 

3.  1996  Distinguished  Faculty  Award  for  Commitment  to  Undergraduate  Research  -  awarded 
by  the  Student  Council  of  the  Johns  Hopkins  University. 

4.  National  Science  Foundation  (NSF)  -  Research  Initiation  Award  -  August  1993.  For 
details  see  Grants  and  Contracts. 

Awards,  fellowships 

Teaching 

1.  Kenan  Award  for  Innovative  Projects  in  Undergraduate  Education  -  May  1995. 
Development  of  an  Experiment  for  Investigation  of  Forced  Convection  Heat  Transfer . 
Awarded  by  the  Kenan  Fund  of  the  Johns  Hopkins  University. 

2.  Subcommittee  for  Electronic  and  Distance  Education  -  SEDE  and  the  Provost’s 
Educational  Technology  Award  Program:  C.  Herman  and  M.  Karweit,  January  1996.  For 
details  see  under  Grants  and  Contracts. 

Research 

1.  NATO  fellowship  for  participation  as  invited  lecturer  in  the  NATO  Advanced  Study 
Institute  on  Low-Temperature  and  Cryogenic  Refrigeration  -  Fundamentals  and  Applications, 
June  23  -  July  5,  2002,  Altin  Yunus-£esme,  Izmir,  Turkey. 

2.  Visiting  Associate  Professor,  Universite  Pierre  et  Marie  Curie  (Paris  6),  Paris,  France,  June 
1- July  3,  2000. 

3.  Chief  of  Naval  Research  Visiting  Scientists  Program:  fellowship  for  the  visit  to  the  Jamie 
L.  Whitten  National  Center  for  Physical  Acoustics,  The  University  of  Mississippi,  Dec.  7-1 1, 
1998. 

4.  NATO  fellowship  for  participation  in  the  NATO  Advanced  Study  Institute  on  Energy 
Conservation  through  Heat  Transfer  Enhancement  of  Heat  Exchangers,  May  25  -  June  5, 
1998,  Cesme/Izmir,  Turkey. 

5.  NATO  fellowship  for  participation  in  the  NATO  Advanced  Study  Institute  on  Cooling  of 
Electronic  Systems,  June  21  -  July  2,  1993,  Cesme/Izmir,  Turkey. 

5.2  Awards:  Martin  Wetzel 

1.  Martin  Wetzel:  Deutscher  Akademischer  Austauschdienst  (DAAD  -  German  Academic 
Exchange  Organization)  fellowship  for  graduate  studies  at  the  Johns  Hopkins  University, 
1994/95  and  1995/96. 

2.  Martin  Wetzel:  ASME  Student  Design  Paper  Award  for  paper  C3, 1996. 

6.  PROJECT  PARTICIPANTS 

6.1  Post-doctoral  scientists 


1 .  Dr.  Yuwen  Chen,  Modeling  of  heat  exchangers  for  thermoacoustic  refrigerators.  Sept.  1998 
-Aug.  1999. 

2.  Dr.  Zdenek  Travnicek,  Design  of  miniature  thermoacoustic  refrigerators,  Nov.  2000  — 
April  2001,  permanent  position  with  the  Institute  of  Thermomechanics,  Academy  of  Sciences 
of  the  Czech  Republic. 

6.2  Graduate  students 

Ph.D,  Students 

1.  Martin  Wetzel,  Thermo-fluid  mechanic  study  of  thermoacoustic  devices,  earned  M.S.  degree 
in  Spring  1997,  completed  Ph.D.  July  1998,  currently  with  BMW  Research,  Munich 
Germany. 

Students  who  have  earned  their  M.S.  degree 

1.  Martin  Wetzel,  Research  topic:  Thermo  fluid  mechanic  study  of  thermoacoustic  devices, 
graduated  in  Spring  1997. 

2.  Ozan  Tutunoglu,  Research  topic :  Design  of  a  miniature  thermoacoustic  refrigerator , 
graduated  Spring  2000. 

Graduate  students  -  in  collaboration  with  European  universities 

(Diploma  Thesis  at  these  universities  approximately  corresponds  to  a  Masters  Thesis) 

I  Technical  University  of  Munich.  Germany: 

1.  Christoph  Bartscher,  Experimental  investigation  of  thermoacoustic  phenomena,  graduated 
June  1994. 

2.  Michal  Volejnik,  Experimental  investigation  of  a  thermoacoustic  refrigerator,  graduated 
February  1995. 

3.  Joachim  Wagner,  Experimental  investigations  of  the  oscillating  flow  and  temperature  fields 
in  a  thermoacoustic  refrigerator,  graduated  January  1 996. 

4.  Markus  Mohne,  Experimental  study  of  heat  transfer  in  a  thermoacoustic  refrigerator ; 
graduated  in  Jan.  1997. 

5.  Andreas  Geishauer,  Experimental  measurements  of  temperature  and  pressure  fields  in  a 
thermoacoustic  refrigerator,  graduated  August  1997. 

13  Fachhochschule  Furtwangen.  Germany: 

1.  Daniel  Herrmann,  Design,  instrumentation  and  initial  experiments  on  a  thermoacoustic 
stack  and  its  heat  exchangers,  graduated  Sept.  1997. 

II  University  of  Lyon.  France: 

1 .  Magali  Brouillat,  Pressure  and  temperature  measurements  in  a  thermoacoustic  refrigerator 
model,  graduated  in  Dec.  Oct  1996. 

7.  ORGANIZATION  OF  CONFERENCES  AND  CONFERENCE 
SESSIONS 

1 .  Session  Organizer 

The  32nd  National  Heat  Transfer  Conference,  Baltimore,  MD,  August  8-12,  1997 

Session  I:  Emerging  Refrigeration  Technologies ,  C.  Herman  Chair,  O.  Knio  Co-Chair. 

Session  H:  Transport  Phenomena  in  Oscillatory  Flows,  C.  Herman  Chair,  O.  Knio  Co-Chair. 

2.  Session  Chair  at  NATO  Advanced  Study  Institute  on  Energy  Conservation  through  Heat 
Transfer  Enhancement  of  Heat  Exchangers,  May  25  -  June  5, 1998,  Cesme/Izmir,  Turkey. 

3.  Session  Organizer  and  Session  Chair 

2  sessions  at  the  137th  Meeting  of  the  Acoustical  Society  of  America,  Berlin,  13-17  March, 
1999 


Session:  Physical  Acoustics:  Thermoacoustics,  6  invited  lectures,  16  papers. 

4.  The  33rd  National  Heat  Transfer  Conference,  Albuquerque,  NM,  Aug.  15-17,  1999 
Session  Organizer:  Session  I:  Emerging  Refrigeration  Technologies 

Session  H:  Transport  Phenomena  in  Oscillatory  Flows 

8.  TECHNOLOGY  TRANSFER 

The  design  optimization  algorithm  developed  by  Wetzel  and  Herman  (1997)  was  used  as  the 
starting  point  for  the  development  of  the  MEMS  thermoacoustic  refrigerator  in  the  project 
sponsored  by  the  Rockwell  Science  Corporation. 


